EFFECTS OF FINITE-PRECISION ARITHMETIC ON 
INTERIOR-POINT METHODS FOR NONLINEAR PROGRAMMING 

STEPHEN J. WRIGHT* 

Abstract. We show that the effects of finite-precision arithmetic in forming and solving the 
hnear system that arises at each iteration of primal-dual interior-point algorithms for nonlinear 
programming are benign, provided that the iterates satisfy centrality and feasibility conditions of 
the type usually associated with path-following methods. When we replace the standard assumption 
that the active constraint gradients are independent by the weaker Mangasarian-Fromovitz constraint 
qualification, rapid convergence usually is attainable, even when cancellation and roundoff errors 
occur during the calculations. In deriving our main results, we prove a key technical result about the 
size of the exact primal-dual step. This result can be used to modify existing analysis of primal-dual 
interior-point methods for convex programming, making it possible to extend the superlinear local 
convergence results to the nonconvex case. 
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1. Introduction. We investigate the effects of finite-precision arithmetic on the 
calculated steps of primal-dual interior-point (PDIF) algorithms for the nonlinear 
programming problem 

(1.1) NLP: min (f>{z) subject to g{z) < 0, 

z 

where : R" ^ R and 5 : R" — > R™ are twice Lipschitz continuously differentiable 
functions. Optimality conditions for this problem can be derived from the Lagrangian 
function C{z, A), which is defined as 

m 

(1.2) £{z, A) = (/.(z) + J2 = H^^) + A^ff(z), 

1=1 

where A € R™ is a vector of Lagrange multipliers. When a constraint qualification 
(discussed below) holds at the point z* , first-order necessary conditions for z* to be 



a solution of (1.1) are that there exists a vector of Lagrange multipliers A* € R™ such 



that the following conditions are satisfied for (z, A) = {z* , A*): 

(1.3) /:,(z,A) = V0(z)+V.g(z)A = O, g{z)<0, A > 0, X^g{z) = 0, 



where 



V<?(z) = [Wgi{z), Wg2{z), Wg^iz)] 



The conditions (1.3) are the well-known Karush-Kuhn- Tucker (KKT) conditions. We 
use S\ to denote the set of vectors A* such that (z*, A*) satisfies ( |l.3| ). The primal-dual 
solution set is defined by 

(L4) S = {z*}xSx. 
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This paper discusses local convergence of PDIP algorithms for (1.1), assuming 
that the algorithm is implemented on a computer that performs calculations accord- 
ing to the standard model of floating-point arithmetic. Because of our focus on local 
convergence properties, we assume throughout that the current iterate (2, A) is close 
enough to the solution set S that superlinear convergence would occur if exact steps 
(uncorrupted by finite precision) were taken. In the interests of generality, we weaken 
an assumption that is often made in the analysis of algorithms for (1.1), namely, 
that the gradients of the active constraints are linearly independent at the solution. 
We replace this linear independence constraint qualification (LICQ) with the weaker 
Mangasarian-Fromovitz constraint qualification (MFCQ) ||l^. MFCQ allows con- 
straint gradients to become dependent at the solution, so that the set S\ of optimal 
Lagrange multipliers is no longer necessarily a singleton, though it remains bounded. 
We continue to assume that a strict complementarity (SC) condition holds, that is, 

(1.5) g^{z*) = =^ A* > 0, for some A* £ Sx. 

In the context of rapidly convergent algorithms, the SC condition makes good sense. 
If SC fails to hold, superlinear convergence of Newton-like algorithms does not occur, 
except for specially modified algorithms such as those that identify the active con- 
straints explicitly (see Monteiro and Wright and El-Bakry, Tapia, and Zhang [^). 

The major conclusion of the paper is that the effects of roundoff errors on the 
rapid local convergence of the algorithm are fairly benign. When a standard second- 
order condition is added to the assumptions already mentioned, the steps produced 
under floating-point arithmetic approach S almost as effectively as do exact steps, 
as long as the distance to the solution set remains significantly greater than the unit 
roundoff u. The latter condition is hardly restrictive, since the data errors made in 
storing the problem in a digital computer mean that the solution set is known only 
to within some multiple of u in any case. 

The conclusions about the effectiveness of the computed steps are not obvious, 
because all three formulations of the linear system that must be solved to compute 
the step at each iteration may become highly ill conditioned near the solution. Our 
analysis would be significantly simpler if we were to make the LICQ assumption 
because, in this case, one formulation of the linear equations remains well conditioned, 
and stability of the three standard formulations can be proved by exploiting their 
relationship to this system of equations. 

This work is related to earlier work of the author on finite-precision analysis of 
interior-point algorithms for linear complementarity problems and linear pro- 
gramming The existence of second-order effects gives the analysis here a 
somewhat different flavor, however. In addition, we go into more depth in checking 
that the computed iterates can continue to satisfy the approximate centrality con- 
ditions usually required in primal-dual algorithms, and in deriving expressions for 
the rate at which the computed iterates approach the solution set. Related work by 
Forsgren, Gill, and Shinnerl |^ deals with one formulation of the step equations for 
the nonlinear programming problem — the so-called augmented form treated here in 
Section ^ — but makes assumptions on the pivot sequence that do not always hold in 
practice. M. H. Wright recently presented an analysis of the condensed form of 
the step equations discussed in Section ^ under the assumption that LICQ holds, and 
found that the computed steps were more accurate than would be expected from a 
naive analysis. 

For linear programming, the PDIP approach has emerged as the most powerful of 
the interior-point approaches. The supporting theory is strong, in terms of global and 
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local convergence analysis and complexity theory (see the bibliography of Wright |26 ). 
Implementations yield better results than pure-primal or barrier-function approaches; 
see Andersen et al. 0. Strong theory is also available for these algorithms when 
applied to convex programming, in which 0(-) and gi{-), i = 1, . . . , m are all convex 
functions; see, for example, Wright and Ralph |3|] and Ralph and Wright [2^ . 
The latter paper drops the LICQ assumption in favor of MFCQ, making the local 
theory stronger in one sense than the corresponding local theory for the sequential 
quadratic programming (SQP) algorithm. The use of MFCQ complicates the analysis 
considerably, however; under LICQ, the implicit function theorem can be used to prove 
a key technical result about the length of the step, while more complicated logic is 
needed to derive this same result under MFCQ. 

A significant by-product of the current paper is to prove the key technical result 
about the length of the rapidly convergent step (Corollary 4.3) under MFCQ and SC, 
even when the problem (|l.l|) i s not convex. This allows the local convergence results 
of Ralph and Wright [|3l|, |2l|, ^ to be extended to general nonconvex nonlinear 
problems. 

The analysis of this paper could also be applied to the recently proposed stabilized 
sequential quadratic programming (sSQP) algorithm (see Wright and Hager |l5| ) , 
in which small penalties on the change in the multiplier estimate A from one iteration 
to the next ensure rapid convergence even when LICQ is relaxed to MFCQ. A finite- 
precision analysis of the sSQP method appears in Section 3.2], but only for the 
augmented form of the step equations. Analysis quite similar to that of the current 
paper could be applied to show that similar conclusions continue to hold when a 
condensed form of the step equations is used instead. We omit the details. 

The remainder of this paper is structured in the following way. Section contains 
notation, together with our basic assumptions about (1.1) and some relevant results 
from the literature. Section |^ discusses the primal-dual interior-point framework, 
defining the general form of eac h ite ration and the step equations that must be solved 
at each iteration. Subsection 3.2 proves an important technical result about the 
relationship between the distance of an interior-point iterate to the solution set S and 
a duality measure /U. Section || describes perturbed variants of the linear systems that 
are solved to obtain PDIP steps, and proves our key results about the effect of the 
perturbations on the accuracy of the steps. 

Section ^ focuses on one form of the PDIP step equations: the most compact 
form in which most of the computational effort goes into factoring a symmetric pos- 
itive definite matrix, usually by a Cholesky procedure. We trace the effect on step 
accuracy of errors in evaluation of the functions, formation of the system, and the 
factorization/solution process. Further, we show the effects of these inaccuracies on 
the distance that we can move along the steps before the interiority condition is vi- 
olated, and on various measures of algorithmic progress. An analogous treatment of 
the augmented form of the step equations appears in Section ^ The conclusions of 
this section depend on the actual algorithm used to solve the augmented system — it is 
not sufficient to assume, as in Section ||, that any backward-stable procedure is used 
to factor the matrix. (We note that similar results hold for the full form of the step 
equations, but we omit the details of this case, which can be found in the technical 
report Q .) We conclude with a numerical illustration of the main results in Section ^ 
and summarize the paper in Section H. 
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2. Notation, Assumptions, and Basic Results. We use B to denote the set 
of active indices at z*, that is, 

(2.1) 6 = {i = l,2,...,m|5,(z*) =0}, 
whereas M denotes its complement 

(2.2) AA={l,2,...,m}\6. 
The set B+ C S is defined as 

(2.3) B+ = {ieB\\*>0 for some A* satisfying (1.3)}. 
The strict complementarity condition ( |l.5| ) is equivalent to 

(2.4) B+ = B. 

We frequently make reference to submatrices and subvectors corresponding to the 
index sets B and A/". For example, the quantities Ag and gei^) are the vectors con- 
taining the components Xi and gi{z), respectively, for i B, while V(7b(z) is the 
matrix whose columns are Vgi{z), i £ B. 

The Mangasarian-Fromovitz constraint qualification (MFCQ) is satisfied at z* if 
there is a vector y £ R" such that 

(2.5) "^gBiz^y < 0. 

The following fundamental result about MFCQ is due to Gau vin |pl| . 

Lemma 2.1. Suppose that the first-order conditions (l.S) are satisfied at z = z* . 
Then S\ is bounded if and only if the MFCQ condition ^2.^ is satisfied at z* . 

This result is crucial because it allows our (local) analysis to place a uniform 
bound on all A in a neighborhood of the dual solution set Sx ■ 

The second-order condition used in most of the remainder of the paper is that 
there is a constant ^ > such that 

(2.6) w'^C,,{z*,X*)w>^\\w\\^ 
for all A* G S\ and all w satisfying 



Vgi{z*)^w = 0, for all i € B+, 
S/g-{z*Yw < 0, for all i G B\B+. 



When the SC condition (L5) (alternatively, (2^)) is satisfied, this direction set is 
simply null V(7|3(z*)"^. 

A simple example that satisfies MFCQ bu t not LICQ at the solution, and that 
satisfies the second-order conditions (^^), (2.7) and the SC condition is as follows: 

(2.8) min zi subject to (zi - 1/3)^ + z^ < 1/9, (zi - 2/3)^ + z^ < 4/9. 

The solution is z* — 0, and the optimal multiplier set is 

(2.9) 5a = {A > 0|2Ai +4A2 = 3}. 

The gradients of the two constraints are the solution are (—2/3,0)^ and (— 4/3, 0)"^, 
respectively. They are linearly dependent, but the MFCQ condition (2.5) can be 
satisfied by choosing y = (1, 0)"^. 
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We use u to denote the unit roundoff, which we define as the smallest number such 
that the following property holds: When x and y are any two floating-point numbers, 
op denotes — , x, or /, and fl{z) denotes the floating-point approximation of a real 
number z, we have 



(2.10) 



fl{x op y) = [x op y){l + e), 



< u. 



Modest multiples of u are denoted by 6^- 

We assume that the problem is scaled so that the values of g and cj) and their first 
and second derivatives in the vicinity of the solution set S, and the values (z, A) them- 
selves, can all be bounded by moderate quantities. When multiplied by u, quantities 
of this type are absorbed into the notation 6^ in the analysis below. 

Order notation O(-) and 8(-) is used as follows: If v (vector or scalar) and e 
(nonnegative scalar) are two quantities that share a dependence on other variables, 
we write v = 0(e) if there is a moderate constant (3i such that ||ti|| < /3ie for all values 
of e that are interesting in the given context. (The "interesting context" frequently 
includes cases in which e is either sufficiently small or sufficiently large, but we often 
use V = 0{fi) to indicate that ||w|| < for all sufficiently small /i that satisfy u, 
for some /3i; see later discussion.) We write w = 8(e) if there are constants (3i and /3o 
such that /3oe < ||u|| < (3ie for all interesting values of e. Similarly, we write v = 0(1) 
if < f3i, and v = 6(1) if /3o < ||w|| < /3i. 

We use the notation S{z,\) to denote the distance from (z, A) to the primal-dual 
solution set, that is, 



(2.11) 



S(z, A) min 



||(z,A)-(z*,A*) 



It is well known (see, for example. Theorem A.l of Wright [g5|]) that this distance can 
be estimated in terms of known quantities at (z. A). We state this result formally as 
follows. 



Theorem 2.2. Suppose that the first- order conditions {l.c ), the MFCQ condition 
{ 2j_i ) and the second-order conditions (2j_6), (2/7) are satisfied at the solution z* . Then 



*/ 2i 0, we have 



(2.12) 



5{z,\) = 6 



Cz{z,\) 
min(A, — g(z)) 



We write the singular value decomposition (SVD) of the matrix V(7g(z*) of first 
partial derivatives as follows: 



(2.13) 



Vgeiz*) = [U V] 





■ 

















where the matrices [ [/ V \ and [ f7 y ] are orthogonal, and S is a diagonal 
matrix with positive diagonal elements. 

Note that the columns of U form a basis for the range space of Vgsiz*), while the 
columns of V form a basis for the null space of V5ib(z*)^. Similarly, the columns of U 
form a basis for the range space of WgB{z*)'^ , while the columns of V form a basis for 
the null space of Vgisiz*). These four subspaces are key to our analysis, particularly 
the subspace spanned by the columns of V. For the computational methods used to 
solve the primal-dual step equations discussed in this paper, the computed step in 
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the yS-components of the muhiphers — that is, AAg — has a larger error in the range 
space of V than in the complementary subspace spanned by the columns of U. The 
errors in the computed primal step Az, the computed step of the A/'-components of 
the multipliers Xj\f, and the computed step in the dual slack variables (defined later) 
are typically also less significant than the error in V'^AXb- We show, however, that 
the potentially large error in AXb does not affect the performance of primal-dual 
algorithms that use these computed steps until /i becomes similar to u^/^. 

When the stronger LICQ condition holds, the matrix V is vacuous, and the SVD 



(2.13) reduces to Vge(z*) — U'SU'^ . Much of the analysis in this paper would simplify 
considerably under LICQ, in part because AAg — the step component with the large 
error — is no longer present. 

We use (Tmin(') to denote the smallest eigenvalue, and cond(-) to denote the con- 
dition number, as measured by the Euclidean norm. 

3. Primal-Dual Interior-Point Methods. 

3.1. Centrality Conditions and Step Equations. Primal-dual interior-point 
methods are constrained, modified Newton methods applied to a particular form of 
the KKT conditions ( |1.3| ). By introducing a vector s £ R™ of slacks for the inequality 
constraint, we can rewrite the nonlinear program as 

min (/)(z) subject to g{z) + s = 0, s > 0, 

{z,s) 



and the KKT conditions (1.3) as 

(3.1) A(2,A)=0, g{z) + s 



0, 



(A,s)>0, X^s = 0. 



Motivated by this form of the conditions, we define the mapping J-{z, A, s) by 



(3.2) 



J^(z, A,s) 



dcf 



V^(z) + Vg(z)A 
g{z) + s 
SAe 



where the diagonal matrices S and A are defined by 

S =^ diag(si, S2, ■ • ■ , Sm), A =^ diag(Ai, A2, . . . , A™), 
and e is defined as 

(3.3) e = (l,l,...,l)^. 

The KKT conditions ( |3.l| ) can now be stated equivalently as 

(3.4) T{z,X,s) = 0, (s,A)>0. 

Primal-dual iterates (z. A, s) invariably satisfy the strict bound (s, A) > 0, while 
they approach satisfaction of the condition J-'{-) = in the limit. An important 
measure of progress is the duality measure /i(A, s), which is defined by 



(3.5) 



/i(A, s) XF s/m. 



When ^ is used without arguments, we assume that /i — /i(A, s), where (z,A, s) is 
the current primal-dual iterate. We emphasize that is a function of (z. A, s), rather 
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than a target value explicitly chosen by the algorithm, as is the case in some of the 
literature. 

A typical step (Az, AA, As) of the primal-dual method satisfies 



(3.6) 



Az 
AA 
As 



-Tiz, A, s) 



where t defines the deviation from a pure Newton step for T (which is also known as 
a "primal-dual affine-scaling" step). The vector t frequently contains a centering term 
tr/xe, where a is a centering parameter in the range [0, 1]. It sometimes also contains 
higher-order information, such as the product AAagAiSaff e, where AAaff and AS'aff 
are the diagonal matrices constructed from the components of the pure Newton step 
(Mehrotra [|9j). In any case, the vector t usually goes to zero rapidly as the iterates 
converge to a solution, so that the steps generated from (3.6) approach pure Newton 
steps, which in turn ensures rapid local convergence. Throughout this paper, we 
assume that t satisfies the estimate 



(3.7) 



All our major results continue to hold, with slight modification, if we replace ( p.7[ ) by 
t — 0{^'^), for some cr e (1,2]. Our essential point remains unchanged; the theoretical 
superlinear convergence rate promised by this choice of t is not seriously compromised 
by roundoff errors as long as fi remains significantly larger than the unit roundoff u. 
To avoid notational clutter, however, we analyze only the case (0). 

Using the definition (|1.2D, we can write the system (|3.6|) explicitly as follows: 



(3.8) 



Czz{z,X) 

V.g(z)^ 




Vff(z) 


s 





/ 

A 



Az ' 




Cz{z, A) 


AA 




g{z) + s 


As 




SAe + t 



Block eliminations can be performed on this system to yield more compact formula- 
tions. By eliminating As, we obtain the augmented system form, which is 



(3.9) 



^zziz,X) ^g{z) 





■ Az ' 






AA 





-Cz{z,X) 
-g{z)+K-H 



By eliminating AA from this system, we obtain a system that is sometimes referred to 
as the condensed form (or in the case of linear programming as the normal equations 
form), which is 



(3.10) 



[Czz{z, A) + Vg{z)kS-^Vg{z)'^] Az 
-Cz{z,X) - Vg{z)KS-^[g{z) - k-H]. 



We consider primal-dual methods in which each iterate (z. A, s) satisfies the fol- 
lowing properties: 



(3.11a) 

(3.11b) 
(3.11c) 



lk/(2;,A)|| < C/i, where r/(z. A) =^ £^(z. A), 

\\rg{z, s)|| < C/i, where rg(z, s) =^ g{z) -f s, 
(A, s) > 0, XiSi > 7/i, for alH = 1, 2, . . . , m. 
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for some constants C > and 7 £ (0,1), where fi is defined as in ( p.5[ ). (In much 
of the succeeding discussion, we omit the arguments from the quantities /i, and 
Tg when they are evaluated at the current iterate (ZjA, s).) These conditions ensure 
that the pairwise products SiXi, i = 1, 2, . . . , m are not too disparate and that the 
first two components of in (3^) can be bounded in terms of the third component. 
They are sometimes called the centrality conditions because they are motivated by 
the notion of a central path and its neighborhoods. Conditions of the type ( 3.11 ) are 
imposed in most path-following interior-point methods for linear programming (see, 
for example, p^). For nonlinear convex programming, examples of methods that 
require these conditions can be found in Ralph and Wright [^l[ |l|, ^ . In nonlinear 
programming, we mention Gould et al. (see Algorithm 4.1 and Figure 5.1) and 
Byrd, Liu, and Nocedal Q. In the latter paper, (3.11a) and (3.111) are imposed 
explicitly, while (3.11c) can be guaranteed by choosing = (1 — 7)/^- Even when the 



choice = /i is made, as in the bulk of the discussion in W, their other conditions 



concerning positivity of (s. A) can be expected to produce iterates that satisfy (3.11c) 
in practice. 

For points (z,A, s) that satisfy ( |3.1lD , we can use fi to estimate the distance 
S{z,X) from (z,X) to the solution set S (see (2.11)). These results, which are proved 
in the following subsection, can be summarized brie fl y as follows. When the MFCQ 
condition (2^) and the second-order conditions (|2.6| ), ( |2.7| ) are satisfied, we have that 
6{z,X) = 0{ij}/'^). When the strict complementarity assumption (L5) is added, we 
obtain the stronger estimate 5(z, A) = 0{ijl). We can use these estimates to obtain 
bounds on the elements of the diagonal matrices S, A, S'~^A, and K~^S in the systems 
above; these bounds are the key to the error analysis of the remainder of the paper. 

3.2. Using the Duality Measure to Estimate Distance to the Solution. 

The main result of this section. Theorem |3.3| , shows that under certain assumptions, 
the distance (5(z, A) of a primal-dual iterate (z. A, s) to the solution set iS can be 
estimated by the duality measure ^. We start with a technical lemma that proves the 
weaker estimate (5(2;, A) = 0(/^^/^). Note that this result does not assume that the 
SC condition ( pTsj ) holds. 

Lemma 3.1. Suppose that z* is a s oluti on of (1.1) at which the MFCQ condition 

(2.1) and the second-order conditions (2^), (12. Ij ) are satisfied. The n for all (z. A) 



with A > for which there is a vector s such that (z,A, s) satisfies ( 3.11 ), we have 
that 



(3.12) 



J(z,A) = 0(//2). 



Proof. We prove the result by showing that 



(3.13) 



C.{z,X) 
min(A, -g{z)) 



and then applying Theorem |2.2| . Since Cz{z,X) = rf = 0{fj,), the first part of the 
vector satisfies the required estimate. For the second part, we have from ( [3.11b ) that 



-ff(^) 



0(m), 



and hence that 



(3.14) min(-5i(z), A^) = min(si, Ai) + 0{fi). 
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Because of (3^) and ( ^.llc ), we have that s^A^ < rrifi and (Aj,Sj) > 0. It follows 
immediately that mm(Xi, Si) < (to/x)^/^ for i — 1,2, ... ,m. Hence, by substitution 
into ( 3.14 ), we obtain 

min(-g,(z). A,) < {mfiY^^ + 0{fi) = 0{fi'^^). 

We conclude that the second part of the vector in ( 3.13| ) is of size 

0(/ii/2), so the 

proof is complete. □ 

The following examples show the upper bound of Lemma is indeed achieved 
and that it is not possible to obtain a lower bound on S{z, A) as a strictly increasing 
nonnegative function of fi. To demonstrate the first claim, consider the problem 

subject to —z < 0. 



1 2 

mm ^z 



The point (z, A, s) — (e, e, e) satisfies 

C,{z,\)^0, g{z) + s = 0. 



so that the conditions ( 3.11 ) are satisfied. Clearly the distance from the point (z. A) 
to the solution set S — (0,0) is \/2e = V2/^^^^. For the second claim, consider any 
nonlinear program suc h th at B = {1,2, . . . , m} (that is, all constraints active) and 
strict complementarity (|l.5| ) holds at some multiplier A* . Then for appropriate choices 
of 7 and C, the point 



(3.15) 



{z,X,.s)^{z*,X*,imfi)/{e^X*)e) 



satisfies the definition ( |3.5| ) and the condition ( 3.11 ) for any fi > 0. On the other 
hand, we have 6{z, A) = S{z* , A*) = by definition, so there are no /3 > and cr > 
that yield a lower bound estimate of the form d{z. A) > /J/i*^. 

We now prove an extension of Lemma 5.1 of Ralph and Wright pl| , dropping the 
monotonicity assumption of this earlier result. 

Lemma 3.2. Suppose that the conditions of Lemma 3.1 hold and in addition that 
the SC condition (1.5) is satisfied. Then for all (z,A, s) satisfying (3.11), we have 
that 



(3.16a) 
(3.16b) 



i(zB^s, = e(fi), A, = 9(1), 
ieJV^s, = 9(1), A, = Q(fi). 



Proof By boundedness of S (Lemma 2.1), we have for all (z,A, s) sufficiently 
close to S that 



(3.17) 



A, = 0(1), s, = -g,{z) + (rg), = 0(1). 



Given (z. A, s) satisfying ( 3.11 ), let P{X) be the projection of A onto the set S\, a nd 
let A* € S\ be some strictly complementary optimal multiplier (for which (L5) is 
satisfied) . From Lemma \A we obtain 



(3.18) 



Using this observation together with smoothness of 0(-) and g(-), we have for the 
gradient of C that 

A(z,A)-£,(z*,A*) 
= V0(z) - V0(z*) + V.g(z)A - V5(z*)A* 

- 0{^^/^) + V.9(z)[A - P(A)] + [Vg(z) - Vg(z*)]P(A) + V.g(z*)[P(A) - A*]. 
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Since P(A) and A* are both in S\, we find from (L3) that the last term vanishes. From 
( 3.18 ) and P(A) — 0(1), the second-to-last term has size 0{^^/^). For the remaining 
term, we have V g{z) = 0(1), and ||A — -P(A)|| < 5{z, A) = 0{fi^^^). By assembling all 
these observations, and using Cz{z* , A*) = 0, we obtain 



(3.19) 



£,(z. A) = A(z, A) - C,{z*, A*) = Oiii'/^). 



Using again that S/g{z*)[P{X) — A*] = 0, we have from ( p.l8[ ) that 

[P(A) - yfW) giz*)] = [PW X*fN9{z*fiz ~ z*) + Oi\\z z*|p) 
(3.20) =0(||z-z*|n =0(Ai). 

By gathering the estimates i ^.l2\ ), ( |3.18| ), ( |3.19D , and ( |3.20D , we obtain 



z 


— z* 


T 


■ C,{z,X)- C,{z\X*) ' 




X 


-A* 








-giz) + giz*) 






z — 


z* 




T 


■ C,{z,X)~C,{z*,X*) 




L A-P(A) J 




-giz) + giz*] 





(3.21) 



+ [P{X}-X*f[-g{z) + g{z*)] 
0(,5(z,A))0(///2) + 0(^)=0(^). 



z — z* 


T 






z — z* 


T r 


A- A* 








A- A* 





By substituting from ( ^.11 ) and using ( ^.21 ), we obtain 

C,iz,X)-C,{z*,X*) 

-giz) + giz*) 

and therefore 

(A - A*)^(s - s*) = -{z - z*)'^rf + (A - X*)'^rg + O(^). 
By using the conditions (3.11a), ( [3.111: ), and the definition (3.5), we obtain 



0(/i), 



A* Si — ^ AjS* 

i=l i=l 

: -(A*)^s - A^s* - -A^s + 0(m) + Oi\\z - z*\\\\rf\\) + 0{\\X - X*\\\\rg\\) = O(^). 



Since (A, s) > and (A*, s*) > 0, all terms X*Si and A^s*, i = 1, 2, . . . , to are nonneg- 
ative, so there is a constant Ci > such that 

< X*Si < Cin, < XiS* < Ci^, for alH = 1, 2, . . . , to. 

For all i €z B, we have A* > by our strictly complementary choice of A*, and so 



(3.22) 



< < ^^i^ < . 

Aj mmigg A, 



Cl dcf „ 

H = C2At. 



On the other hand, we have by boundedness of S\ and our assumption (3.11c) that 



(3.23) 



Si> — > 7minM, for all 1 = 1,2,, 
A, 



,m, 
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for some constant 7i„in > 0. By combining ( 3.22 ) and ( 3.23 ), we have that 

Si — 0(/i), for all i E B. 
For the Ag component, we have that 

SiXi>'-ff^ Xi> — >-^, for alH e 
Si 62 



Hence, by combining with (3.17), we obtain that 

A, = 6(1), for aU i e B. 



This completes the proof of (3.16a). We omit the proof of (3.16b), which is similar. □ 
Next, we show that when the strict complementarity assumption is added to the 
assumptions of Lemma 3.1, the upper bound on the distance to the solution set in 
(3.12) can actually be improved to 0(/x). 



Theorem 3.3. Suppose that z* is a solution of (1.1) at which the MFCQ con- 
dition (2.5), the second-order conditions ^.(\ ), (2.'/), and the SC condition (|77^ are 
satisfied. Then for all (z. A) with A > for which there is a vector s such that (z. A, s) 
satisfies (3.11), we have that 



(3.24) 



<S(z,A)-0(^). 



Proof From (3.11a), we have directly that r/ — 0{fx). We have from (3.11) and 



( p.l6aD that 

5i(z) = -Si + (rg), = 0(//), A, = 6(1), A, > for alH e 

so that 

(3.25) m\n{—gi{z),\i) — —gi{z) — 0{y), for alH e S, 
whenever is sufficiently small. For the remaining components, we have 

(3.26) min(-.g,(z), Ai) = A, = 0(m), for aU i e TV. 



By substituting ( ^.llaD , ( |3.25| ), an d (|3 26| ) into ( |2.12| ), we obtain the resuh. □ 

Similar conclusions to Lemma 3.2 and Theorem 3.3 can be reached in the case 



of linear programming algorithms. The second-order conditions (2.6), (2.7) are not 
relevant for this class of problems, and the SC assumption ( |l.5| ) holds for every linear 
program that has a solution. 

4. Accuracy of PDIP Steps: General Results. By partitioning the con- 
straint index set {1,2,..., m} into active indices B and inactive indices A/", we can 
express the system ( ^.9[ ) without loss of generality as follows: 





>Czz(2;, A) 


^geiz) ygj\f{z) 




Az 




-C,(z,\) 

-gB{z) + A^^tB 


(4.1) 


V.98(^)^ 


-Db 




AAb 






_ Vg^/{zY 


-Du 




AAa/- 




. -gj^{z) + h~^tj^ _ 



where Db and Dj^ are positive diagonal matrices defined by 
(4.2) Db = Agi^B, Dm = A^^Sa^. 
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When the SC condition (L5) is satisfied, we have from Lemma that the diagonals 
of_Dg have size 6(/u) while those of Dj^ have size By elimin ating AA^ from 

(4.1), we obtain the following intermediate stage between (3^) and ( 3.10| ): 



(4.3) 



H{z,X) Vgeiz) 



Az 



-9b{z) + A_g tB 



where we have defined 



(4.4) 



H{z, A) 1^*' CUz, A) + ygu{z)DjlygM{zf 



In this section, we sta rt by proving a key result about the solutions of perturbed 
forms of the system (O). Subsequently, we use this res ult as the foundation for 
proving results about the three alternative formulations (|3.8|), ( |3.9|), a nd ( 3.10 ) of 
the PDIP step equations. The principal reason for our focus on (O) is that the 



proof of the main result can be derived from fairly standard linear algebra arguments. 
Gould ||l^, Section 6] obtains a system similar to (O) for the Newton equations for 
the primal log-barrier function, and notes that the matrix approaches a nonsingular 
limit when certain optimality conditions, including LICQ, are satisfied. Because we 
replace LICQ by MFCQ, the matrix in (4.3) may approach a singular limit in our 
case. 



We note that the form (4.3) is also relevant to the stabilized sequential quadratic 
programming (sSQP) method of Wright |^ and Hager jl^; that is, slight modifica- 
tions to the results of this paper can be used to show that the condensed and aug- 
mented formulations of the step equations for the sSQP algorithm yield good steps 
even in the presence of roundoff errors and cancellation. We omit further details in 
this paper. 

Errors in the step equations arise from cancellation and roundoff errors in evalu- 
ating both the matrix and right-hand side and from roundoff errors that arise in the 
factorization/solution process. We discuss these sources of error further and quantify 
them in the next section. In this section, we consider the following perturbed version 
of (gj): 



(4.5) 



H{z,\) + En VgB{z) + Ei2 
S/gB{zV + E2i -DB + E22 



w 

y 



ri 



VgB{z*yr3+ri 



Here, E is the perturbation matrix (appropriately partitioned and not assumed to 
be symmetric) and ri, r^, and represent components of a general right-hand side. 
Note the partitioning of the second right-hand side component into a component 
WgB{z*)^r3 in the range space of Vgg(z*)-^ and a remainder term 7-4. When LICQ is 
satisfied, the range space of V(?e(z*)-^ spans the full space, so we can choose r4 to be 
zero. Under MFCQ, however, we have in general that must be nonzero. The main 
interest of the results below is in isolating the component of the solution of (L5) that 
is sensitive to r^. 

To make the results applicable to a wider class of linear systems, we do not 
impose the assumptions that were needed in the preceding section to ensure that the 
matrices Db and defined by (4.5) have diagonals of the appropriate size. Instead, 
we assume that the diagonals have the given size, and derive the application to the 
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linear systems of interest (those that arise in primal-dual interior-point methods) as 
a special case. 



Our results in this and later sections make extensive use of the SVD ( 2.13 ) of 
Vgg(2:*). They also make assumptions about the size of the smallest singular value 
of this matrix, and about the size of the smallest eigenvalue of V'^£zz{z* , X*)V , the 
two-sided projection of the Lagrangian Hessian onto the active constraint manifold. 

Theorem 4.1. Let (z, A) be an approximate primal-dual solution of (l.l) with 
S{z, A) — 0{fi), and suppose the diagonal matrices Djg and Dj^ defined by (^.i ) have 



all their diagonal elements of size 6(/i). Suppose that the perturbation submatrices in 
( R^ satisfy 



(4.6) 



Ell = Su/ti + 0{p), E2i,Ei2, E22 



and that the following conditions hold for some f3 > 

(4.7a) 
(4.7b) 

(4.7c) cr 



u/^<l, U<1, 
o-min(S) > /3max(^^/^,u/^), 
{V'^Cz,{z*,X*)V) > /3max(/^i/3,u/^), for aU A* e Sx- 



Then if (3 is sufficiently large (in a sense to be specified in the proof), the step {w,y) 
computed from (■(■•-) satisfies 



iU^y,V^w,U^w) = 0{\\ri\\ + \\rs\\ + \\r4), 

F^y = 0(||ri|| + ||r3|| + ||r4||^). 



Proof. If we define 

yu = u^y, 



yv 



we have 

y = Uyu + Vyv, 

Using this notation, we can rewrite (E^) as 



w = U Wjj + VWy. 



(4.8) 



U^MuU U^MiiV U^Mi2U U^Mi2V 

V^MiiU V^MiiV V^Mi2U V^Mi2V 

U^M2iU U'^M2iV U^M22U U^M22V 

V^M2iU V'^M2iV V^M22U V'^M22V 









Wy 




yu 




yv 



where we have defined 



(4.9) 



Mil = H{z, A) + Ell, M12 = V.9b(z) + E12, 
M2i^VgB{zf + E21, M22 = -Db + E22, 



and H{-, ■) is defined in (4^). From ( 2.13 ), we have 

V^VgB{z*f^G, U^VgB{z*f 
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The fact that V'^ annihilates '^gei^*)'^ is crucial, because it causes the term with 
to disappear from the last component of the right-hand side of ( [l.§| ), which becomes 

From the definitions (HJ) and Q), the perturbation bound ( |4.6| ), our assumptions 
that Djj^ = 0{fi) and S{z,X) = 0{fj,), compactness of S, and the fact that Czz is 
Lipschitz continuous, we have that 



(4.10) 



(4.11) Mn = Czz{z\X*) + 5 J 11 + 0{fi), 

for some A* E S\. Using these same facts, we have likewise that 



so by substituting from (2.13), we have that 

(4.12a) U^M2iU = T: + du + 0{fi), U'^ M21V = 6^ + 0{fi), 

(4.12b) V'^M^iU = Su + 0{fi), V'^A'hiV ^S^ + 0{ii). 

Similarly, from the definition of A/12, we have 

(4.13a) U'^MuU = J: + S^ + 0{fi), U'^ M12V = 6^ + 0{fi), 

(4.13b) V'^Mi2U = Su + 0{fi), V'^Mi2V = S^ + 0{n). 

For the A/22 block, we have from ( p^ ) and (4.6) that 

(4.14a) U'^M22U = -U^DbU + 5u = 0(m) + '^u, 

(4.14b) U^M22V = 0{ii) + (5u, V'^M22U = 0{ii) + S^, 

(4.14c) V^M22V = -V^DbV + Su^ Mvv + Su, 



where Mvv = -V^DbV has all its singular values of size 0(/i), so that 



(4.15) 



vv 



e(/.-i) 



Using these estimates together with ( 4.10| ), we can rewrite ( |4.8| ) as 
(4.16) 



Q 

Mvv 



Ell E12 
E21 E22 







Wy 




yu 




yv 





V^Ti 



where 
(4.17) 



Q 



(4.18) 



dcf 



E 
d^u/M + o(m) 5^1 ^i + 0{p) 

<5u + o{^l) 

Nuu Nuv Si 
Nyu Nvv 
E2 
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while 
(4.19) 

and 
(4.20) 



Ell 



El2, E21 





+ 0(m) 

5u + 0{fi) (5u + 0(m) 



E. 



22 



For purposes of specifying the required conditions on /? in ( 1.7b ) and ( 4.7c ), we 
define k to be a constant such that expressions of size 5^^ and 0(/i) that arise in the 
perturbation terms in the coefficient matrix in (4.16) can be bounded by ku and K/i, 
respectively. For example, we suppose that the perturbations in Ei, E2, and Nyv 
can bounded as follows: 



(4.21a) 
(4.21b) 

and that 



ISl - Ell < K(/i + u), \\T.2-^\<K{^1 + M) 

\\Nvv - V^CU^\X*)V\\ < Kiu/fx + y), 



(4.22)||^ii|| < k(u + ^), 11^1211 < k(u + /x), ||i;2i|| < k(u + ^), 11^2211 <'«u. 

From ( [4.21aD and ( H.7b| ), we have that 

||Si-S|| < Km&^{^i^'\u| ^x) < (K//3)a,„i„(S) < (k//3)||E||. 

It is therefore easy to show that if f3 can be chosen large enough that /3 > 2k (while 
still satisfying (4.7b) and (4.7c)), then 



(4.23) ||Ei||<2||E| 
Similarly, we can show that 

(4.24) ||E2||<2||E| 



|Sr'll<2||I]-i| 



E^^ll < 2\\T.- 



(4.25) ||7Vyy|| < 2||f^/:,,(z*,A*)l>||, IliVv^^ll < 2\\{V^C,,{z*,X*)V)-'\\. 

Note, too, that because of Lipschitz continuity of Czz and compactness of S, and the 
bounds (1.7a), the norms of Njjij, Njjv, Nvu, Nvv, and E are all 0(1). Hence the 
matrix Q is itself invertible, and we have 



(4.26) = 

Noting that 
(4.27) 












vv 



-l^i^NuvNy^ -ti\Nuu - NuvNylNvu)^2^ 



(Q + Ell 



(I + Q-'Eiir^Q~\ 



we examine the size of Q ^Eii- Note first from ( 4.7b| ) and ( 4.7c| ) together with ( 1.23| ), 
( p4| ), and ( p^ ) that 



(4.28a) 
(4.28b) 



|E7^II< 4(u/m)-\ ||E2-1|| <|(u^)-\ IITV^^II < |(u^)-\ 



/3 



f3' 



2 



1/3 



|E-i|| < 



AT-l II < ±,,-1/3 
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By forming the product of ( 4.26| ) with ( 4.19 ) and using the bounds in ( 4.28 ), we ca n 
show that the norm of Q^^En can be made less than 1/2 provided that (3 in ( 4.71 ), 
( 4.7c ) is sufficiently large. The (3,3) block of Q~^Eii, for instance, has the form 

-Sr'iVc/yiV^7^((5u + 0(a*)) + ^i\NiJu - NuvN-l,Nvu)^2'iS^ + 0{p)). 



Because of (4.22), its norm can be bounded by a quantity of the form 

Cn{\\t^'\\ WN-l^W + pr^ll WNylW + Pr^ll IIS^-^II) ((u/^)/i + ^) , 



(for some C that depends on ||£22(2*, A*)||), which in turn because of ( [4.28[ ) is bounded 
by the following quantity: 



Provided that f3 is large enough that this and the other blocks of Q^^E u can be 
bounded appropriately, we have that < 1/2, and therefore from ( 4.27 ) we 

have 

||(Q + ^n)-i|| =2||Q-i|!. 

Our conclusion is that for /? satisfying the conditions outlined in this paragraph, the 
inverse of the (1, 1) block of the matrix in ( 4.1(j) can be bounded in terms of 
which because of ( |42^ ), ( ^1 , ( ji^ ), and (H.26| ) can in turn be bounded by a finite 
quantity that depends only on the problem data and not on fi. 
Returning to ( 4.16| ), and using ( 4.20| ), we have that 



Wy 

yu 



(4.29) 



-(Q + En)-'Ei2yv + {Q + Eii)- 

0(|l^i2|l|lyy||) + 0(|lri|l + |lr3|l + ||r4| 
0(Ai)||yy||+0(||ri|| + ||r3|| + ||r4||). 



Meanwhile, for the second block row of (4.16), we obtain 



(4.30) 



yv 



-{Myv + E22) ^E2l 



yu 



+ {Mvv + E22)''V' r4 



Since from ( |4.15| ), ( ^.20D , and ( [4.7a| ), we have 

{Mvv + E22)-' = {1 + MylE22)-'Mvl = (/ + SJ^i)Myl, = Oip.'^), 



it follows from ( ^.30D and ( [4.20| ) that 
yv = 0{^i-')0{pi) 

By substituting from ( [l.29| ), we obtain 

\\yv\\ ^ 0{fi)\\yv\\ + 0{\\n\ 



yu 



+ 0it,'')0{\\r4 



+ 11^311 + ||r4||) + 0(||r4||/Ai), 
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and therefore 



\yv\ 



Odlnl 



'■3 



as claimed. The estimate for {wjj,Wy,yu) is obtained by substituting into (4.29). □ 



The conditions (4.7) need a httle explanation. For the typical value u = 10^^^^, the 
minimum value of the quantity max(/j,-'^/^, u//Lt) is 10^"', achieved at /i^^^. Moreover, 
we have max(/i^/'^, u/^) < 10~^ only for fx in the range [10"^"*, 10~^1. It would seem, 



then, that the problem would need to be quite well conditioned for (4.71:) and (4.7c) 
to hold and that n may have to become quite small before the results apply. We 
note, however, that the purpose of the bounds ( 4.7b| ) and (4.7c) is to ensure that 
the inverse of Q + En can be bounded independently of ji, and that for this purpose 
they are quite conservative. That is, we would expect to find that ||(Q + is 
not too much larger than the norm of the inverse of the corresponding exact matrix 
(the first term on the right-hand side of ( 4.17 )) for ^ not much less than the smallest 
eigenvalues of S and Cy_z{z* ,X*)V . 

The requirement that u//i and ^ both be small in ( [4.7| ) may not seem to sit well 
with expressions such as 0(/x) and 0(/x^), which crop up repeatedly in the analysis 
and which assert that certain bounds hold "for all sufficiently small /i." As noted in 
the preceding paragraph, this requirement implies that the analysis holds for in a 
certain range, or "window," of values. Similar windows are used in the analysis of 
S. Wright [g^ 11^, and M. H. Wright and numerical experience indicates that 
such a window does indeed exist in most practical cases. We expect the same to be 
true of the problem and algorithms discussed in this paper. 

At this point, we assemble the assumptions that are made in the remainder of 
the paper into a single catch-all assumption. 

Assumption 4.1. 

(a) z* is a solution of (1-1), so that the condition (l.S) holds. The MFCQ con- 



dition {2.L), the second-order conditions (2.6), (2.7), and the SC condition 
(l.L) are satisfied at this solution. The current iterate (z,\s) of the PDIP 



algorithm satisfies the conditions (3.11), and the right-hand side modification 
t satisfies (3.1). 

The quantities fi, u (2.1L ), Czz{z* , A*), Yi, and V (2.1c) satisfy the conditions 



(b) 

(kji- 

From our observations following 



(4.31) 



Dr- 



we have under this assumption that 



Our next result considers a perturbed form of the system (4.1), with a general 
right-hand side. By eliminating one component to obtain the form (4.3), we can apply 
Theorem 4.1 to obtain estimates of the dependence of the solution on the right-hand 



side components. 
Theorem 4.2. 



Suppose that Assumption ^.1 holds. Consider the linear system 





- Ell 


VgB{z) + Ei2 


VgAr{z) + Ei^ ' 




w 




V- E21 


—Djs + E22 


E23 




y 


. ^9n{zV - 


f 


E32 


-D^f + Es:^ 




q 



(4.32) 



?'5 

ygi3{z*)^re 
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where 

(4.33a) 
(4.33b) 



Eii=5u/H, E33 = Su/fJ-'^, 
El2; E21, E22 — Su, -E13, £"31, i?23i £-32 



(5u//i- 



Then the step (w,y,q) satisfies the following estimates: 

(C/^y,u;) = 0(||r5|| + ||r6|| + ||r7||+M||r8||), 

V^y = Odirsll + llrell + ||r7||^ + + 0(M))||r8||), 

q = 0(m) [\\r4 + \\r,\\ + \\rs\\] + i6J„ + 0{„mrj\\. 



that 



Proof. From (4.31) and the assumed bound (4.33a) on the size of i?33, we have 



(4.34) = -(/ - D]^'E33)-'D]^' = {I + 0{pi)S^/pi^)0{pi) = 0{fi). 

By eliminating q from ( 4.32| ), we obtain the reduced system 



H{z,X) + En VgB{z) + Ei2 
\79b{zV + E2i -D13 + E2 



^22 



w 

y 



rr, + Oi^l)\\rs\\ 
S/gB(z*)^re+rT + SJrs\\ 



where from ( |4.7D and (4.4), we obtain 

En = En - {\/g^f{z) + El3){-D^f + E33)-\^g^f{z f + E31) - ^gMz)D]^^^g^f{zf 
= (5u/m + 0(/x), 

^12 = E12 - (V5A^(z) + Ei3)hD^ + E33)-^E32 = 5^ + 0(1)0(m)<5u/m = ^u, 
E21 = E21 - E23{-Dm + E^zV^-^guizf + £;3i) = <5u + {5^/ ^l)0{^i)0{l) = 5u, 

E22 = E22 - E23{~D^ + E33y^E32 = 5u + {SJfifOifl) = Su- 



These perturbation matrices satisfy the assumptions of Theorem O, which can be 
applied to give 

(4.35a) {U^y, V^w, U^w) = 0(||r5|| + llrgH + ||r7|| + ^i\\rs\\), 

(4.35b) y^y^Odlrsll + ||r6|| + ||r7||//i) + (<5u//i + 0(^))||r8||. 

From the last block row of and using ( |l3^ ), ( ^35| ), we obtain 

q = (-Dm + E33y^ [rs - {Vg^izf + E3i)w ~ E32y\ 
= 0{p) [||r8|| + ||u;|| + (5u^)||2/||] 
^0{^i) [||r5|l + ||r6|| + ||r7|| + ||r8||] + 

'^u [llrsll + llrell + ||r7||//i + (5u//i + 0(M))lk8ll] 
= 0{^i) [llrsll + llrell + \\rs\\] + {SJ^i + 0(A^))||r7|| . 

□ 

An estimate for the solution of the exact system (3^) follows almost immediately 
from this result. This is the key technical result used by Ralph and Wright [^l|, ^ to 
prove superlinear convergence of PDIF algorithms for convex programming problems. 
The result below, however, does not require a convexity assumption. 
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Corollary 4.3. Suppose that Assumption ^^(a) holds. Then the (exact) solu- 
tion (Az, AA, As) of the system j^S.^ satisfies 

(4.36) (Az, AA, As) = O(^). 



Proof. Note first that Assumption 4.1(b) holds trivially in this case for n suffi- 
ciently small, because our assumption of exact computations is equivalent to setting 
u = 0. We prove the result by identifying the system (|4.l|) with (4.32) and then 



applying Theorem 4.2 



For the right-hand side, we note first that, beca use of smoothness of g, Taylor's 
theorem, the definition (2T) of S, and Theorem |3.3| , 

geiz) - gsiz*) + VgB{z*f{z - z*) + 0{\\z - z*f) 
(4.37) ^VgB{z*f{z^z*) + 0{^l''). 



We now identify the right-hand side of (4T) with ( 4.32 ) by setting 

rb = -Cz{z,X), 
re = {z- z*), 

r? = -\/gB{z*f{z - z*) - gB{z) + K^Hb, 
rs = -9J\r{z) + Kjftu. 

The s izes of these vectors can be es timated by using ( 3.11 ), Lemma 3.2, ( 4.37 ), The- 
orem 3.3, and the assumption (3.7) on the size of t to obtain 

(4.38) r-5 = 0(M), r6 = 0(M), = O(m'), rg = 0(1). 

(By choosing re and r^ in this way, we ensure that the terms involving ||r7||//x in the 
estimates of the solution components in Theorem 4^ are not grossly larger tha n the 
other terms in these expressions.) We complete the identification of (4T) with ( 4.32| ) 
by setting all the perturbation matrices En, E12, • ■ • , £"33 to zero and by identifying 
the solution vector components Az, AAg, and AXj\/ with w, y, an d q, respectively. By 
directly applying Theorem L2, substituting the estimates ( [4.3^ ), and setting 5u = 



(since we are assuming exact computations), we have that 

{U^AXb, Az) = 0(/i), ^^^AAb = 0{^l), AX^f 



o(m). 



To show that the remaining s olut ion component As of (3.8) is also of size 0{pL), 

as 



we write the second block row in (3 

As = -{g{z) + s) 



Vg{z)^Az, 



from which the desired estimate follows immediately by substituting from ( 3.11b ) and 
Az = 0{n). □ 



The next result uses Theorem [4.2| to compare perturbed and exact solutions of 
the system of the system (4T). 

Corollary 4.4. Suppose that Assumption ^.1 holds. Let (w,y,q) be obtained 
from the following perturbed version of ( jj'.4 )-' 



Czz{z,X) - 

^9b{zV- 



En V<?b(2) 

£^21 
- £31 



E 



12 



E 



13 



(4.39) 



V5a/'(2;) 
Db + £22 E23 
E32 —Dj\f + £33 

-A(z,AH/i 
-gBiz) + As^tB + f2 
-9N-{z) + A^^A/- -f /a 



w 

y 
q 
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where Eij, i,j — 1,2,3, satisfy the conditions (4-.3c ) and fi, /2, and are all of size 
Su- Then i/(Az,AA, As) is the (exact) solution of the system we have 

{Az~w,U^iAXB-y))^S^, 

V^{AXB-y) = 5Jfi, 
AX^r - q = Su- 



Proof. By combining (4.39) with (fill), we obtain 





^zz(z, X) - 




^98{z) - 


f E^2 ^9m{z) - 


\- Ei3 




w- A 


z 




Vg6(z)^H 








E22 


E23 






V - AXb 




. ^gM^V - 






E32 






E33 




_ q - AA^ _ 








' h ' 




Ell E12 


E13 




Az 




(4.40) 








h 




E21 E22 


E23 




AXb 










h 




E31 E32 


E33 




AXat 





From the bounds on the perturbations E in (4.33) and the resuh of Corollary 4.3 
have for the right-hand side of this expression that 



we 





r5 


dot 


" /l " 




Ell E12 E13 




Az 










/2 




E21 E22 E23 




AXb 












. /3 . 




E31 E32 E33 














S 


u + ((^u/m)/-* + (^uM + (f^u/^)/^ 






(4.41) 








Su- 


l-(5uM + ^uM+ {Su/l^)lJ. 










_ 


- {Su/fJ.)fJ. + (<5u/m)m + (^u/m^)^^ . 







Using these estimates, we can simply apply Theorem O to ( 4.40| ) (with re = 0) to 

obtain the result. □ 

For later reference, we show how the estimates of Corollary 4.4 can be modified 
when the perturbations have a special form. Suppose that 

(4.42) E23^0, S33-'^u/m, f2 = Uf^ + 0{ii^), where /f = 5u, 

where U is the matrix from ( ^.13| ). Instead of setting rg as in the proof above, 
we set 

(using ( ^.13 ) to obtain an rg for which Vgg(z*)-^r6 — U f^)- B y mo difyin g ( 4.41 ) to 
account for the remaining perturbations, we can identify (4.40) with (4.32) by setting 



dcf 



(4.43) 



/l 




Ell E12 E13 




A 


z 


/2 - Uf^ 




E21 E22 E23 




AXb 


f3 




E31 E32 E33 




_ AX^ 


Su + (<5u/M)Ai 


+ + (<5u/Ai)M 












Su + {5u/ti)lJ- + 


{5^/^l)^i + {5^/^J^)^l _ 







O(m') 

Su/lJ- 



Using these modified right-hand side estimates, we can apply Theorem 4.2 to obtain 
the following improved bound on one of the components: 



(4.44) 



V'^iAXB-y)^Oi^l). 
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The bounds on the other components remain unchanged. 

We emphasize that the conditions ( [3.11 ), and in particular ( 3.11c ), are critical 
to the results of this and all the remaining sections of the paper. These conditions 
enable Lemma 3.2, which in turn enable us to assert that the diagonals of Djs all have 
size 6(/i) while those of Dj^ all have size Q{^~^) (see ( 4.3l| )). This neat classification 
of the diagonals of D into t wo ca tegories drives all the subsequent analysis. The 
motivation for conditions like (3.11) in the literature for path- following methods (with 
exact steps) is not unrelated: It allows us to obtain bounds on the steps and to show 
that we can move a significant distance along this direction while ensuring that (3.11) 
continues to be satisfied at the new iterate. (See, for example, pq, Chapters 5 and 
6] and its bibliography for the case of linear programming and | |31| , for the 

case of nonlinear convex programming.) In the analysis above, we obtain bounds on 
the errors (rather than the steps themselves) when perturbation terms of a certain 
structure appear in the matrix and right-hand side. 

Many practical implementations of path-following methods for linear program- 
ming do not explicitly check that the condition (3.11c) is satisfied by the calculated 
iterates (see, for example, Q and Q). However, the heuristics for "stepping back" 
from the boundary of the nonnegative orthant by a small but significant quantity are 
motivated by this condition, and it is observed to hold in practice on all but the most 
recalcitrant problems. 

5. The Condensed System. Here we consider an algorithm in which the con- 
densed linear system ( ^.10 ) is formed and solved to obtain Az, and the remaining 
step components A A and As are recovered from (3.8). We obtain expressions for the 



errors in the calculated step (Az, AA, As) and discuss the effects of these errors on 
certain measures of step quality. We also derive conditions under which the Cholesky 
factorization applied to ( |3.1C| ) is guaranteed to run to completion. 
Formally, the complete procedure can be described as follows: 

procedure condensed 

given the current iterate (z. A, s) 

form t he co efficient matrix and right-hand side for ( |3.10 ); 
solve (|3.10|) using a backward stable algorithm to obtain Az; 
set AA = D-^[g{z) - A-H + Vg(z)'^Az]; 
set As = -{g{z) + s) - Vg(z)^Az. 



We have used the definition 
the system ( |3.10 ) here as follows: 



of the matrix D. For convenience, we restate 



(5.1) [£,,{z,X) + Vg{z)D-^Vg{zf] Az = -C,{z,X) - S/g{z)D-^[g{z) - A~h]. 



Note that this procedure requires evaluation of _D ^ — S ^A, rather than D itself. 
5.1. Quantifying the Errors. When implemented in finite-precision arith- 



metic, solution of (5.1) gives rise to errors of three types: 

- cancellation in evaluation of the matrix and right-hand side; 

- roundoff errors in evaluation of the matrix and right-hand side; 

- roundoff errors that accumulate during the process of factoring the matrix 
and using triangular substitutions to obtain the solution. 

Cancellation may be an issue in the evaluation of the nonlinear functions Czz{z, A), 
Cz{z,X)^ g{z), and \7g{z), because intermediate terms computed during the additive 
evaluation of these quantities may exceed the size of the final result (see Golub and 
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Van Loan |12, p. 61]). The intermediate terms generally contain rounding error (which 
occurs whenever real numbers are represented in finite precision). Cancellation be- 
comes a significant phenomenon whenever we take a difference of two nearly equal 
quantities, since the error in the computed result due to roundoff in the two argu- 
ments may be large relative to the size of the result. If, as we can reasonably assume, 
intermediate quantities in the calculations of our right-hand sides remain bounded, 
the absolute size of the errors in the result is 6^- In the case of Cz{z,X) and geiz), 
the final result in exact arithmetic has size 0{fi), so that the error takes on a large 
relative significance for small values of /i. This fact causes the error bound in some 
components of the solution to be larger than in others, as we see in ( |3.6c ) below. In 
summary, the computed versions of the quantities discussed above differ from their 
exact values in the following way: 



(5.2a) 


computed Czziz, X) ^ C 


,z{z,X) + F, 


(5.2b) 


computed Cziz, A) ^ Cz{z, X) + f 


(5.2c) 


computed ^g{z) ^ ^g{z) + F = 


■ Vgf3(z) " 
. V5A/-(z) _ 


+ 


Fb 

F^ \ ' 


(5. 2d) 


computed g{z) <~ g{z) + / 


9b{z) 

_ miz) 


+ 





where F, /, F, and / are all of size ^u- Earlier discussion of cancellation in similar 
contexts can be found in the papers of S. Wright ||, H, |^ and M. H. Wright [||. 
The second source of error is evaluation of the matrix D^^. From the model 



(2.10) of floating-point arithmetic and the estimates (3.16) of Lemma 3.2, we have 
that 



(5.3a) 
(5.3b) 



computed 1?^^ ^ {Db + Gb) ^, Gg = ^(5u, 
computed DJ^ ^ {D^f + Ga/-)^\ G^s — S^/ fJ., 



where Gb and Gj\f are both diagonal matrices that can be composed into a single 
diagonal matrix G. 

Third, we account for the error in forming the matrix and right-hand side of 



(5.1) from the computed quantities described in the last two paragraphs. Since we 
are now dealing with floating-point numbers, the model ( |2.10|) applies; that is, any 
additional errors that arise during the combination of these floating-point quantities 
have size u relative to the size of the result of the calculation. Since the norm of 
the coefficient matrix is of size 0{fi~^) while the right-hand side has size 0(1) (see 



(3.11)), we represent these errors by a matrix F of size Su/fJ- and a vector / of size 
Finally, we account for the error that arises in the application of a backward-stable 



method to solve (5.1). Specifically, we assume that the method yields a computed 
solution that is the exact solution of a nearby problem whose data contains relative 
perturbations of size u. The absolute sizes of these terms would therefore be (5u//i in 
the case of the matrix and in the case of the right-hand side. Since these errors 
are the same size as those discussed in the preceding paragraph, we incorporate them 
into the matrix F and the vector /. 



Summarizing, we find that the computed solution Az of (5.1) satisfies the follow- 
ing system: 



(5.4) 



Czz{z, X)+F+ (V.g(z) + F)iD + G)-\Vg(z) + Ff + F Az 
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-A(z, A) - / - {Vg{z) + F){D + G)-'[g(z) + f- A-H] + /, 



where the perturbation terms F, F, F, G, /, /, and / are described in the paragraphs 
ab ove. By " unfo lding" this system and using the partitioning of F, G, and / defined 
in ( |3.2[ ) and (^.3| ) , we find that Az also satisfies the following system, for some vectors 
y and q: 



(5.5) 



Vg^izf + Fi 





+ Fb VgM{z) - 






' Az ' 




Gb 






y 







-Dm- 


Gm 










-£,{z,X)-f- 










-9b{z) + A^^tB 


-fB 















This system has precisely the form of ( 4.39| ) (in particular, the perturbation matrices 
satisfy the appropriate bounds). Hence, by a direct application of Corollary 4.4, we 
conclude that 



(5.6a) 
(5.6b) 
(5.6c) 



Az — Az = (5u, 
U^{A\b - 2/) - <5u, 
V^{A\B-y)^6Jti. 



We return now to the recovery of the remaining solution components AA and As 
in the procedure condensed. We have from Assumption 4.1 together with (3.111: ), 
Lemma O, (Ell), (|^, (|53a|), (O), and (Op that 



{5.7a.)gB{z) = rg{z,s)B-SB = 0{fi), Ag^ = 9(1), Az = Az + = 0{fi), 
(5.7b) {Db + Gb)-' = {I + D^'Gb)D^' = {I + S^)-^0{^l-') = OUr'). 

Since t ~ 0(/i^), we have from our model ( ^.10 ) that the floating-point version of the 
calculation of AAg in the procedure condensed satisfies the following: 



AAg = {Db + Gb)-' gB{z) + /k - Ag^g + (V.gf3(z) + Fb)^Az + fiS^ 



+ 5u- 



(The final term i5u arises from ( 2.10| ) because our best estimate of the quantity in the 
brackets at this point of the analysis is 0(/i), so from (5.7b) the result has size 0(1).) 
Meanwhile, we have from the second block row of (5.£) that 



y={DB + Gb)-' gB{z) + Ib - ^b'^b + {^9b{z) + Fb)'^Az 



By a direct comparison of these two expressions, and using {Db + Gb) ' — 0{ijl), we 
find that 



(5.8) 



AAg-y 



By combining (5.8) with (5.6b) and (5.6c), we find that 

(5.9) U^{A\b - A\b) - V^{A\b - AAg) = 5^|^Ji. 



24 



STEPHEN J. WRIGHT 



F or the "nonb asic" part AAtv", we have from ( ^.llbD , Lemma U ( p6|) , ( F6a| , 
( F3b|) , and ( ^^Sll) that 

(5.10a) gM{z)^0{l), ^J}^0{^^^), K^z = 0{ii), 

(5.10b) {Dj^ + G^)-i = (/ + DJ^Gj^r^Djl = DJ^' + fiS^ = 0{fi). 

By using tj\f = 0(/i^) and applying the model ( 2.1C| ) to the appropriate step in the 
procedure condensed, we obtain 



AAaa = [Dm + Gj^)-^ \gM{z) + fM- ^J-tM + (VgA^(z) + Fu fK'z + 5^] + 
By comparing this expression with the corresponding exact formula, which is 



and by using the bounds (5.1C) and the fact that and have size ^u, we obtain 

AX^f - AX^f = ^iSu + {Dj^ + Gj^y^[fM + 6^] + 

[{Dj^ + Gm)-^ - DJ^%_^{z) - A^^W] + 
{D^f + G^)-\VgM(z) + F^fAz - Dj^VgHizfAz 
= f^S^ + {D;^ + GMT^'^gAzfi^z - Az) + ^^^A^] 
[{D^f + Gu)-^ - Dj;/]Vg^{zfAz 
(5.11) =Ai<^u. 

Finally, for the recovered step As, we have from the last step of procedure con- 
densed, together with ( |3.11l| ), ( |5.2dD , ( |5.7bD , and ( ^.lOD that 



As 



-{g{z) + f + s)- (V.g(z) + i^)^Az + 5^, 



where the final term accounts for the rounding error ( |2.1(]| ) that arises from accumu- 
lating the terms in the sum, which are all bou nded. By s ubstit uting the expression 
for the exact As together with the estimates ( 5. 2d ) and ( 5.71 ) on the sizes of the 
perturbation terms, we obtain 

As = -(.g(z) + s) - V5(z)'^Az - / - V.g(z)^(A^ - Az) - F'^Az + ii5^ 
(5.12) =As + (5u. 

We summarize the results obtained so far in the following theorem. 

Theorem 5.1. Suppose that Assumption ^.1 holds. Then when the step 
(Az, AA, As) is calculated in a finite- precision environment by using the procedure 
condensed ( and where, in particular, a backward stable method is used to solve the 
linear system for the Az component) , we have that 



(5.13a) 
(5.13b) 
(5.13c) 



(Az - Az, C/^(AAb - AAh), As - As) = 5^, 

V'^ {AXrs - AXb) = 5J ii, 
AXj^f - AAa/- = ^^u. 
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This theorem extends the result of M. H. Wright |23 for accuracy of the computed 
sohition of the condensed system by relaxing the LICQ assumption to MFCQ. When 
LICQ holds, the matrix V is vacuous, so the absolute error in all components is of size 
at most Su- The higher accuracy ( 5.13c ) of the components AAa/" (also noted in ||2^ ) 
does not contribute significantly to the progress that can be made along the inexact 



direction (Az, AA, As), in the sense of Section 5.3. 

We return briefly to the case discussed immediately after Corollary 



4.4, in which 



the perturbations have the special form ( 4.42 ), using these results to show that the 
bound (5.13b) can be strengthened when /e satisfies 



(5.14) 



This case is of interest when the ca ncella tion errors in computing geiz) are smaller 
than the estimate we made following ( 5.2(3 ), possibly because of use of higher-precision 
arithmetic or the fact that the computation did not require differencing of quantities 
whose size is large relative to the final result. When ( 5.14 ) holds, we see by comparing 
( p9|) with dJ) that 



£^23 = 0, E33=G^ = SJfx, f2 = UU^fB + 0{p^), where /e = 5u. 



Therefore, we deduce from ( 4.44|) that ( 5.6c ) can be replaced by 

^^(AAb -y) = 0(/i). 
Using ( ^.8[ ) and /i ^ (5u, we can therefore replace ( 5.13b| ) in this case by 
(5.15) V^{AXb - AXb) ^ 0{n). 

5.2. Termination of the Cholesky Algorithm. In deriving the estimate 



( p.6[ ), we have assumed that a backward stable algorithm is used to solve (5.1). Be- 
cause of (p.(]|), (p.7D, and the SC condition, and the estimates of the sizes of the 



diag onals of D (from (4.2) and Lemma 3.2), it is easy to show that the matrix in 
( ^.l[ ) is positive definite for all sufficiently small fi. The Cholesky algorithm is there- 
fore an obvious candidate for solving this system. However, the condition number of 
the matrix in (5^) usually approaches oo as /i | 0, raising the possibility that the 
Cholesky algorithm may break down when fj, is small. A simple argument, which we 
now sketch, suffices to show that successful completion of the Cholesky algorithm can 
be expected under the assumptions we have used in our analysis so far. 

We state first the following technical result. Since it is similar to one proved by 
Debreu Theorem 3], its proof is omitted. 

Lemma 5.2. Suppose that M and A are two matrices with the properties that M 
is symmetric and 



Mx > allMll 



for some constant a > 0. Then for all /i such that 

_dcf . (a\\A\\^ \\A\\ 
< ^ < ^ = mm f 



we have that 



i\\M\\ ' a||M||y ' 



x'^{M + fi-^AA^)x > -\\xf, for aU x. 
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We apply this result to by setting 

M = C,,{z, A) + WgM{z)Dj/Wgj^{zf = C,,{z, A) + 0(/i), 



-1/2 



(where again we use (4.2) and Lemma 3.2 to derive the order estimates). The con- 
ditions (2^), (2/7), and strict compl eme ntarity ensure that this choice of M and A 



satisfies the assumptions of Lemma 5.2. The result then implies that the smallest 
singular value of the matrix in (5.1) is positive and of size 6(1) for all values of ^ 
below a threshold that is also of size 0(1). Since D — 0(/i^^), the largest eigenvalue 
of this matrix is of size 0{^,^^), so we have 



(5.16) 



cond{C^^{z,\) + V g{z)D-^V g{z)^) = 0(p-^). 



(An estimate similar to this is derived by M. H. Wright |23, Theorem 3.2], under the 
LICQ assumption.) It is known from a result of Wilkinson (cited by Golub and Van 
Loan jl2| p. 147]) that the Cholesky algorithm runs to completion if (7„(5uCond(-) < 1, 
where g„ is a modest quantity that depends polynomially on the dimension n of the 
matrix. By combining this result with (5.16), we conclude that for the matrix in (pj). 



we can expect completion of the Cholesky algorithm whenever /i 3> ^u- That is, no 
new assumptions need to be added to those made in deriving the results of earlier 
sections. 

We note that this situation differs a little from the case of linear programming 
where, because second-order conditions are not applicable, it is usually necessary to 



modify the Cholesky procedure to ensure that it runs to completion (see [30 



5.3. Local Convergence Virith Computed Steps. We begin this section by 
showing how the quantities rf, Vg, and /i change along the computed step (Az, AA, As) 
obtained from the finite-precision implementation of the procedure condensed. We 
compare these with the changes that can be expected along the exact direction 
{Az, AA, As). We then consider the effects of these perturbations on an algorithm of 
the type in which the iterates are expected to satisfy the conditions (3.11). Rapidly 



convergent variants of these algorithms for linear programming problems usually al- 
low the values of C and 7 in these conditions to be relaxed, so that a near-unit step 
can be taken. We address the following question: If similar relaxations are allowed in 
an algorithm for nonlinear programming, are near-unit steps still possible when the 
steps contain perturbations of the type considered above? 

We show in particular that for the computed search direction, the maximum step 
length that can be taken without violating the nonnegativity conditions on A and s 
satisfies 

(5.17) 1 - Amax = (^u/m + 0(Ai), 

while the reductions in pairwise products, fJ-, rf, and Vg satisfy 

(5.18a) (Ai -I- aAXi){si + aAsj) = (1 - a)X,Si + 6^ + 0{n'^), i = 1, 2, . . . , to, 
(5.18b) /i(A + aAX, s + aAs) = (1 - a)fi + 6u + 0{p^), 
(5.18c) r/(z -I- aAz, A + aAA) = (1 - a)rf +5^ + 0(/i^), 
(5.18d) rg{z + aAz, s + aAs) = (1 - a)rg + + 0(^^). 
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The corresponding maximum steplength for the exact direction satisfies 
(5.19) 1 - a,„ax = 0(Ai), 



while the reductions in Vf, Vg, and ^ satisfy 



(5.20a) 
(5.20b) 
(5.20c) 
(5.20d) 



(A, + aAX,){s, + aAs,) = (1 - a)\,s, + 0{p^), 
fi{X + aAX, s + aAs) = (1 - a)n + 0{fi^), 
rf{z + qAz, A + aAX) = (1 - a)rf + 0{i?), 
rg{z + aAz, s + aAs) = (1 - a)rg + 0{i?). 



z = 1,2,. 



, ,m, 



Our proof of the estimates ( 5.17 ) and ( 5.18 ) is tedious but not completely straight- 
forward, and we have included it in the Appendix. 

It is clear from ( 5.17 ) and ( 5.18 ) that the direction (Az,AA, As) makes good 
progress toward the solution set S. If the actual steplength a is close to its maximum 
value cimaxj in the sense that 

(5.21) Amax - a = (Ju/a* + 0(/^), 
we have by direct substitution in (5.17) and ( ^.18 ) that 

/i(A + aAA, s + aAs) =5^ + 0{fi'^), 
rf{z + aAz, A + aAA) = + O(m^), 
rg{z + aAz, s + aAs) = 6u + 0{n'^). 

These formulae suggest that finite precision does not have an observable effect on the 
quadratic convergence rate of the underlying algorithm until /i drops below about 
y/u. Stopping criteria for interior-point methods usually include a condition such as 
yU < lO^u ov n < ^Jvl (see, for example, so that is not allowed to become so 
small that the assumption /i 3> u made in (1.7) is violated. 

In making this back-of-the-envelope assessm ent, h owever, we have not taken into 
account the approximate centrality conditions (3.11), which must continue to hold 
(possibly in a relaxed form) at the new iterate. These conditions play a central role 
both in the analysis above and in the convergence analysis of the underlying "exact" 
algorithms, and also appear to be important in practice. Typically (see, for example, 
Ralph and Wright [^), the conditions (3.11) are relaxed by allowing a modest increase 
in C and a modest decrease in 7 on the rapidly convergent steps. We show in the next 
result that enforcement of these relaxed conditions is not inconsistent with taking a 
step length a that is close to dmax, so that rapid convergence can still be observed 
even in the presence of finite-precision effects. 

Theorem 5.3. Suppose that Assumption ^.1 holds. Then when the step 
{Az, AX, As) is calculated in a finite-precision environment by using the procedure 
condensed, there is a constant C such that for all r e [0, 1/2] and all a satisfying 

(5.22) ae [0,1-C't-1(u/m + m)], 

the following relaxed form of the approximate centrality conditions holds: 



(5.23a) 
(5.23b) 
(5.23c) 



r/(z + aAz, A + aAA) < C(l + r)//(A + aAA, s + aAs), 
rg{z + aAz, s + aAs) < C(l + r)//(A + aAA, s + aAs), 



{Xi + aAXi){s^ + aAs,) > 7(1 



-r)M(A- 
for all i 



aAX, s + aAs) 
-- 1,2, ...,m. 
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where C is the constant from conditions (3.11). Moreover, when we seta to its upper 
bound in (5.2i ), we find that 



(5.24) 



6{z + aAz, A + aAA) < t^^S^ + 0{f/)). 



Proof. From (3.11) and (5.16), we have that 

\\rf{z + aAz, A + q;AA)|| 
^{l-a)\\rf\\+S^ + 0{fi^) 
< C{l~a)fi + S^ + 0{n^) 

= C(l + t)(1 - a)fi - Ct(1 - a)/i + + 0{i?) 

= C(l + r)/i(A + aAA, s + aAs) - Cr(l - a)^i + 5^ + 0(/i^). 



We deduce that the required condition ( |5.23aD will hold provided that 

5u + 0{^?) < CT{l~a)fi. 

Since by definition we have that 6^ + 0{(j?) < C{u + /i^) for some positive constant 
C, we find that a sufficient condition for the required inequality is that 

(l-a)> (C/C)r-i(u^ + /i), 

which is equivalent to ( 5.22| ) for an obvious definition of C. Identical logic can be 
applied to ||rg|l to yield a similar condition on a. 

For the condition ( ^.23c ), we have from (3.11) and (5.18) that 

{X., + aA\.,)is, + aAs,) 
= (1 - a)X^s,+d^ + 0{fi^) 

> (l-a)7/i + (5u + 0(M') 

= 7(1 - t)(1 - a)^i + 7r(l - a)fi + 6^ + 0{p^) 

= 7(1 - T)^l{\ + aAA, s + aKs) + 7r(l - a)^i + 6^ + 0{fM^). 



Hence, the condition (5.23c) holds provided that 

7T(l-a)M + 5u + 0(Ai^) >0. 

Similar logic can be applied to this inequality to derive a bound of the type ( ^.22 ), 
after a possible adjustment of C. 

Finally, we obtain ( 5.24 ) by substituting a 1 — Ct^^{u/ fi + /i) into ( 5.18 ) and 
apply ing Theorem |]^. (Note tha t, despite the relaxation of the centrality conditions 
( ^.23 ), the res ult of Theor em |3.3| still holds; we sim ply mo dify the proof to replace C 
by (3/2)C in ( |3.11a| ) and ( |3.11b|) , and 7 by 7/2 in ( |3'.llcD .) □ 

6. The Augmented System. In this section, we consider the case in which the 
augmented system (3.9) (equivalently, (4.1)) is solved to obtain (Az, AA), while the 
remaining step component As is recovered from (3.8). The formal specification for 
this procedure is as follows: 
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procedure augmented 

given the current iterate {z, A, s) 



form the coefficient matrix and right-hand side for (4.1); 
solve (^31) to obtain (Az, AA); 
set As = -ig{z) + s) - Vgiz^Az. 



Much of our work in analyzing the augmented system form ( 11 ) ha s already been 
performed in Section ^; the main error result is simply Corollary iA. However, we 
can apply this result only if the floating-point errors made in evaluating and solving 
this system satisfy the assumptions of this corollar y. In particular, we need to show 



that the perturbation matrices Eij, i,j = 1, 2, 3 in ( 4.3£ ) satisfy the estimates ( |4.33| ). 

This task is not completely straightforward. Unlike the condensed and full-system 
cases, it is not simply a matter of assuming that a backward-stable algorithm has been 



used to solve the system (4.1). The reason is that the largest terms in the coefficient 
matrix in ( |3.9| ) — the diagonal elements in the matrix Dj^ — have size 0(/i~^). The 
usual analysis of backward-stable algorithms represents the floating-point errors as a 
perturbation of the entire coefficient matrix whose size is bounded by times the 



matrix norm — in this case, ^. However, Corollary 4.4 requires some elements of the 
perturbation matrix to be smaller than this estimate; in particular, the submatrices 
£'i2, E21, and E22 need to be of size 5^. Therefore, we need to look closely at the 
particular algorithms used to solve (^]^) to see whether they satisfy the following 
condition. 

Condition 6.1. The solution obtained by applying the algorithm in question to 



the system (i.l) in floating-point arithmetic is the exact solution of a perturbed system 



in which the perturbations of the coefficient matrix satisfy the estimates {(.So), while 
the right-hand side is unperturbed. 

We focus on diagonal pivoting methods, which take a symmetric matrix T and 
produce a factorization of the form 



(6.1) PTP^ ^ LYL"^, 

where P is a permutation matrix, L is unit lower triangular, and Y is block diagonal, 
with a combination of 1 x 1 and symmetric 2x2 blocks. The best-known methods of 
this class are due to Bunch and Parlett ||^ and Bunch and Kaufman [|| , while Duff et 
al. and Fourer and Mehrotra |]l^ have described sparse variants. These algorithms 
differ in their selection criteria for the 1x1 and 2x2 pivot blocks. In our case, the 
presence of the diagonal elements of size 0(^~^) (from the submatrix D_\f = Aj^S_\f) 
and their place in these pivot blocks are crucial to the result. 

We start by stating a general result of Higham ||l^ concerning backward stability 
that applies to all diagonal pivoting schemes. We then examine the Bunch-Kaufman 
scheme, showing that the large diagonals appear only as 1 x 1 pivots and that this 
algorithm satisfies Condition |6.l[ (In |l^. Theorem 4.2], Higham actually proves that 
the Bunch-Kaufman scheme is backward stable in the normwise sense, but this result 
is not applicable to our context, for the reasons mentioned above.) 

Next, we briefly examine the Bunch-Parlett method, showing that it starts out 
by selecting all the large diagonal elements in turn as 1 x 1 pivots, before going on to 
factor the remaini ng m atrix, whose elements are all 0(1) in size. This method also 



satisfies Condition 6T . We then examine the sparse diagonal pivoting approach es of 
Duff et al. [0 and Fourer and Mehrotra [|o|, which may not satisfy Condition 6T 



because of the possible presence of 2 x 2 pivots in which one of the diagonals has size 
0(/i~^). These algorithms can be modified in simple ways to overcome this difficulty. 
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possibly at the expense of higher density in the L factor. We then mention Gaussian 
elimination with pivoting and refer to previous results in the literature to show that 



this approach satisfies Condition |6.1| . Finally, we state a result like Theorem |5.3| 
about convergence of a finite-precision implementation of an algorithm based on the 
augmented system form. 

Higham Theorem 4.1] proves the following result. 

Theorem 6.1. Let T be an n symmetric matrix, and let x be the computed 
solution to the linear system Tx — b produced by a method that yields a factorization 



of the form (6.1), with any diagonal pivoting strategy. Assume that, during recovery 
of the solution, the subsystems that involve the 2x2 diagonal blocks are solved via 
Gaussian elimination with partial pivoting. Then we have that 

(6.2) (r + AT)x = 5, \m<S^m + P'^\L\\Y\\L'^\P) + 5l, 

where L and Y are the computed factors, and \A\ denotes the matrix formed from A 
by replacing all its elements by their absolute values. 

In Higham's result, the coefficient of u in the Su term is actually a linear poly- 
nomial in the dimension of the system. The partial pivoting strategy for the 2x2 
systems can actually be replaced by any method for which the computed solution 
of Ry — d satisfies {R + AR)y = d, where R is the 2x2 matrix in question and 
|A_R| < (5u|i?|. This property was also key in an earlier paper of S. Wright p^, who 



derived a result similar to Theorem 6.1 in the context of the augmented systems that 
arise from interior-point methods for linear programming. 

All the procedures below have the property that the growth in the maximum 
element size in the remaining submatrix is bounded by a modest quantity at each in- 
dividual step of the factorization. (In the case of Bunch-Kaufman and Bunch-Par lett, 
this bound averages about 2.6 per elimination step; see Golub and Van Loan Sec- 
tion 4.4.4].) As with Gaussian elimination with partial pivoting, exponential element 



growth is possible, so that L and Y in (6.1) contain much larger elements than the 



original matrix T. Such behavior is, however, quite rare and is confined to pathologi- 
cal cases and certain special problem classes. In our analysis below, we make the safe 
assumption that catastrophic growth of this kind does not occur. 

6.1. The Bunch-Kaufman Procedure. At each iteration, the Bunch- 
Kaufman procedure chooses either a 1 x 1 or 2 x 2 pivot by examining at most two 
columns of the remaining matrix, that is, the part of the matrix that remains to be 
factored at this stage of the process. It makes use of quantities Xi defined by 

Xi = max \Tij\, 

j I 

where in this case T denotes the remaining matrix. We define the pivot selection 
strategy for the first step of the factorization process. The entire algorithm is obtained 
by applying this procedure recursively to the remaining submatrix. 

set 1/ = (1 + %/T7)/8; 

calculate xii store the index r for which xi = 

if |Tii| >z^Xi 

choose Til as a 1 X 1 pivot; 

else 

calculate Xr', 
i{Xr\Tii\>iyxl 
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choose Til as a 1 X 1 pivot; 
else if jr^rl > i^Xr 

choose Trr as a 1 X 1 pivot; 

else 

choose a 2 X 2 pivot with diagonals Tn and Trr', 
end if 
end if. 

For each choice of pivot, the permutation matrix Pi is chosen so that the desired 
1 X 1 or 2 X 2 pivot is in the upper left of the matrix PiTP^ . If one writes 



PiTP^ 



R 

C f 



where R is the chosen pivot, the first step of the factorization yields 



(6.3) 



/ 




■ R 




■ / R-^C^ ' 






T 




I 



where T^T~ CR-^C'^ is the matrix remaining after this stage of the factorization. 

At the first step of the factorization, the quantities xi and Xr (if calculated) both 
have size 0(1), since the large elements of this matrix occur only on the diagonal. 
Since a 2 x 2 pivot is chosen only if 



\Tii\<i^Xi and \Trr\ < vXr, 

it follows immediately that both diagonals in a 2 x 2 pivot must be 0(1) 
pivot chosen by this procedure is one of three types: 



Hence, the 



(6.4a) 
(6.4b) 
(6.4c) 



1x1 pivot of size 0(1); 

2x2 pivot in which both diagonals have size 0(1); 
1x1 pivot of size 8(/i~^). 



In fact, the pivots are one of the types (6.4) at all stages of the factorization, not 
just the first stage. The reason is that the updated matrix T in ( |6.3[ ) has the same 
essential form as the original matrix T — its elements are all of size 0(1) except for 
some large diagonal elements of size Q{^~^). We demonstrate this claim by showing 
that the update CR^^C^ that is applied to the remaining matrix in ( |6.3| ) is a matrix 
whose elements are of size at most 0(1), regardless of the type of pivot, so that it 
does not disturb the essential structure of the remaining matrix. When the pivots are 
of type (6.4a) and ( |6.4b| ), the standard argument of Bunch and Kaufman can be 
applied to show that the norm of CR~^C^ is at most a modest multiple of ||C||. We 
know that ||C|| = 0(1), since C consists only of off-diagonal elements, so we conclude 
that \\CR-^C^\\ = 0(1) in this case as claimed. For the other pivot type ( |6.4c) ), we 
have R = 8(/^^^) and C — 0(1), so the elements of CR~^C^ have size 0{fi), and 
the claim holds in this case too. 



In the rest of this subsection, we show by using Theorem 3.1 that Condition 6.1 
holds for the Bun ch-Kaufman al gor ithm. In fact, we prove a stronger result: When 
T in Theorem 3T is the matrix (4T), the perturbation matrix AT contains elements 
of size 6u, except in those diagonal locations corresponding to th e el ements of D^, 
where they may be as large as Su/fJ-- Given the bound on |AT| in (|6.2|), we need only 



32 



STEPHEN J. WRIGHT 



to show that P-^^lLj |y| |I/|-^P has the desired structure. In fact, it suffices to show 
that the exact factor product P-^|L| |y| |L|-^P has the structure in question, since the 
difference between these two products is just in size. 

We demonstrate this claim inductively, using a refined version of the arguments 
from Higham |]l7| . Section 4.3] for some key points, and omitting some details. For 
simplicity, and without loss of generality, we assume that P = I. 

When n = 1 (that is, T is 1 x 1), we have that L = 1 and Y = T, and the result 
holds trivially. When fi = 2, there are three cases to consider. If the matrix contains 
no elements of size 8(/i^^), then the analysis for general matrices can be used to 
show that = 0(1), as required. If either or both diagonal elements have 

size 9(/i^^), then both pivots are 1x1, and the factors have the form 



(6.5) 



L 



1 



T21/T1 



Y 



Til 






Til/Ti 



Two cases arise 
(i) 



(ii) 



A diagonal of size 0(1) is accepted as the first pivot and moved (if necessary) 
to the (1, 1) position. We then have 

iT^iil > i^Xi = 1^X2 ^ '^\T2i\, 

and therefore iTai/Tnl < l/i^ and hence iTi^/Tn] < \T21\lv = 0(1). If the 
(2, 2) diagonal is also 0(1), we have that L = 0{1) and Y — 0(1), and we are 
done. Otherwise, T22 = 0(^^^), and so the (2,2) element of Y satisfies this 
same estimate. We conclude from (|6.5| ) that |L||F||L|-^ also has an 9(^~^) 
element in the (2, 2) position and 0(1) elements elsewhere. 
A diagonal of size 9(/Lt~^) is accepted as the first pivot and moved (if neces- 
sary) to the (1, 1) position. We then have 

|T2i/rn| =0(/i), \Tl/Tll\^0{^,). 



It follows from (6^) that 



ILiiyiiLi 



\Tii\ 

\T2l\ 



\T2l\ 

\T22\+0{y) 



which obviously has the desired structure. 
We now assume that our claim holds for some dimension n > 2, and we prove 
that it continues to hold for dimension n + 1. Using the notation of (6.3) (assuming 
that Pi = /), and denoting the factorization of the Schur complement T in ( |6.3| ) by 
f = ZFL^, we have that 



(6.6) 



T = LYL' = 



I 




■ R 




■ / R-^C^ 


CR-^ L 




Y 







It follows that 

(6.7) \L\\Y\\L\^ 



\R\ 
\CR-^ 



R\ ICR 



\R\\R-'C^\ 
^^\\R\\R-^C^\ + 



LWYULl 



Since, as we mentioned above, the norm of CR-^C^ is at most 0(1), the Schur 
complement T = f - CR^^C'^ has size 0(1) except for large 8(/i ^) elements in 
the same locations as in the original matrix. Hence, by our inductive hypothesis. 
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|L||y||L|-^ has a similar struc ture , and we need to show only that the effects of the 
first step of the factorization ( |6.3| ) do not disturb the desired structure. 

For the case in which i? is a pivot of type either ( 6.4a ) and ( |6.4b| ), Higham ||l7| , 
Section 4.3] shows all elements ofboth \CR-^\\R\ and \CR-^\\R\\R-'^C'^\ are bounded 
by a modest multiple of either xi (if Tu was selected as the pivot because it passed 
the test |Tii| > i^Xi) or (Xi +Xi'), where r is the "other" column considered during 
the selection process. In our case, this observation implies that both |Ci?^^||i?| and 
|Ci?^^||i?||i?^^C^| have size 0(1). By combining these obse rvati ons wi th th ose of the 
preceding paragraph, we conclude that for pivots of types ( |6.4a ) and ( |6.4b ), "large" 
elements of the matrix in (6.7) occur only in the diagonal locations originally occupied 
by D^f. 

For the remaining case — pivots of type ( 3.4c ) — we have that C has size 0(1) while 
has size 0(u). Therefore, \CR~^\\R\ has size 0(1) and \CR'^\\R\\R-~^C^\ has 



size 0(/i), while \R\, which occupies the (1, 1) position in the matrix (6.7), just as it 
did in the original matrix T, has size 0(/i^^). We conclude that the desired structure 
holds in this case as well. 



We conclude from this discussion that Condition 6.1 holds for the Bunch-Kaufman 



procedure. We show later that the perturbations arising from other sources, namely, 
roundoff and cancellation in the evaluation of the matrix and right-hand side, also 



satisfy the conditions of Corollary 4.4, so this result can be used to bound the error 
in the computed steps. 



Finally, we note that it is quite possible for pivots of types (6.4a) and (3.4b) 
to be chosen while diagonal elements of size 0(/z~^) still remain in the submatrix. 
Therefore, a key assumption of the analysis of Forsgren, Gill, and Shinnerl |, Theo- 
rem 4.4] — namely, that all the diagonals of size <d{^~^) are chosen as 1 x 1 pivots before 
any of the other diagonals are chosen — may not be satisfied by the Bunch-Kaufman 
procedure. 

6.2. The Bunch-Parlett Procedure. The Bunch-Parlett procedure is con- 
ceptually simpler but more expensive to implement than Bunch-Kaufman, since it 
requires 0{n^) (rather than 0(n)) comparisons at each step of the factorization. The 
pivot selection strategy is as follows. 



set iy= (l + V17)/8; 

calculate Xoff = \Trs \ = max^^j |Ty |, Xdiag = l^^ppl = max^ \Tii\; 

if Xdiag > fXoS 

choose Tpp as the 1x1 pivot; 

else 

choose the 2x2 pivot whose off-diagonal element is Trs', 
end if. 



The elimination procedure then follows as in (3.3) 



It is easy to show that the Bunch-Parlett procedure starts by selecting all the 
diagonals of size Q{fi^^) in turn as 1 x 1 pivots. (Because of this property, it satisfies 
the key assumption of mentioned at the end of the preceding section.) The update 
CR~^C^ generated by each of these pivot steps has size only 0{ii), so the matrix 
that remains after this phase of the factorization contains on ly (1) elements. The 
remaining pivots are then a combination of types ( 3.4a ) and ( |6.4b| ). 

By using the arguments of t he preceding subsection in a slightly simplified form, 
we can show that Condition 3.1 holds for this procedure as well. 
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6.3. Sparse Diagonal Pivoting. For large instances of (pTl|), the Bunch- 
Kaufman and Bunch-Parlett procedures are usually inefficient because they do not 
try to maintain sparsity in the lower triangular factor L. Sparse variants of these 
algorithms, such as those of Duff et al. and Fourer and Mehrotra ||l^, use pivot se- 
lection strategies that combine stability considerations with Markowitz-like estimates 
of the amount of fill-in that a candidate pivot will cause in the remaining matrix. 

At each stage of the factorization, both algorithms examine a roster of possible 
1x1 and 2x2 pivots, starting with those that would create the least fill-in. As soon 
as a pivot is found that meets the stability criteria described below, it is accepted. 
Both algorithms prefer to use 1x1 pivots where possible. 

For candidate 1x1 pivots. Duff et al. |fl p. 190] use the following stability criterion: 



(6.8) 



\R- 



where the notation R and C is from (6^) and p G [2,oo) is some user-selected pa- 
rameter that represents the tolerable growth factor at each stage of the factorization. 
For a 2 X 2 pivot, the criterion is 



(6.9) 



\R- 







< 


p 




. llC-,2||oo . 




p 



where C.^i and are the two columns of C. The stability criteria of Fourer and 

Mehrotra |]l^ are similar. 

As they stand, the stability tests (6.8) and ( |6.9| ) do not necessarily restrict the 
choice of pivots to the three types ( |6.4|) . If a 1 x 1 pivot of size 0(/i^^) is ever 
considered for structural reasons, it will pass the test ( |6.8| ) (the left-hand side of this 
expression will have size 0{p)) and therefore will be accepted as a pivot. However, it 
is possibl e th at 2x2 pivots in which one or both diagonals have size Q{f j,~^) may pass 
the test (3^) and may therefore be accepted. Although the test (|6.9| ) ensures that 
the size of the update CR^^C'^ is modest (so that the update T = T — CR~^C^ does 
not disturb the large-diagonal structure of T), there is no obvious assurance that the 
matrix |L||y||L|"^ in (6.7) mirrors the structure of |T|, in terms of having the large 
diagonal elements in the same locations. The terms |Ci?~^||i?| and \CR''^\\R\\R^^C'^\ 



in ( |6.7| ) may not have size 0(1), as they do for pivots of the three types (6.4) arising 
from the Bunch-Kaufman and Bunch-Parlett selection procedures. 

The Fourer-Mehrotra algorithm does, however, rule out the possibility of a 2 x 2 
pivot in which both diagonals are of size Q{fj,~^). It considers a 2 x 2 candidate only if 
one of its diagonal elements has previously been considered as a 1 x 1 pivot but failed 
the stability test. However, if either of the diagonals had been subjected to the test 
( |6.8[ ), they would have been accepted, as noted in the preceding paragraph, so this 
situation cannot occur. 

If the sparse algorithms are modified to ensure that all pivots have one of the 

or (|6.9D, then 



three types (6.4), and all continue to satisfy the stability tests (6.8 
simple arguments (simpler than those of Section |6.1| !) can be applied to show that 
Condition 3.1 is satisfied. One possible modification that achieves the desired efffect 



is to require that a 2 x 2 pivot be allowed only if both its diagonals have previously 
been considered as 1 x 1 pivots but failed the stability test (|6.8D. 



6.4. Gaussian Elimination. Another possibility for solving the system (4A) is 
to ignore its symmetry and apply a Gaussian elimination algorithm, with row and/or 
column pivoting to preserve sparsity and prevent excessive element growth. Such 



FINITE-PRECISION EFFECTS IN NONLINEAR PROGRAMMING 



35 



a strategy satisfies Condition 6A. In |Q, the author uses a result of Higham iQ 
to show that the effects of the large diagonal elements are essentially confined to the 
columns in which they appear. Assuming that the pivot sequence is chosen to prevent 
excessive element growth in the remaining matrix, and using the notation of ( 4.32| ) 
and ( 4.33 ), we can account for the effects of roundoff error in Gaussian elimination 
with perturbations in the coefficient matrix that satisfy the following estimates: 

Ell, E21, E^i, Ei2, E22, = <5u, Ei3, E23, E32 = Su/fJ'- 



These certainly satisfy the conditions ( 4.33 ), so Condition 6A holds. 

6.5. Local Convergence with the Computed Steps. We can now state a 
formal result to show that when the evaluation errors are taken into account as well as 
the roundoff errors from the factorization/solve procedure discussed above, the accu- 
racies of the computed steps obtained from the procedure augmented, implemented 
in finite precision, satisfy the same estimates as for the corresponding steps obtained 



from the procedure condensed. The result is analogous to Theorem 5.1 



Theorem 6.2. Suppose that Assumption ^.1 holds. Then when the step 
(Az, AA, As) is calculated in a finite-precision environment by using the procedure 
augmented, where the algorithm used to solve ^4-^) satisfies Condition 6.1, we have 



(6.10a) 
(6.10b) 
(6.10c) 



(Az - Az, U'^{AXb - AAe), As - As) = S^, 

V^{A\b - AXb) = 6u/ti, 
AA^ — AA^ = Su- 



Pro of. The proof follows from Corollary 44 when we show that the perturbations 
to ( 4T ) from all sources — evaluation of the matrix and right-hand side as well as the 
factorization/solution pro cedu re — satisfy the bounds required by this earlier result. 

Because of Condition 6.1, per turbations arising from the factorization/solution 
procedure satisfy the bounds ( 4.33| ) . The expressions (|]^) show that the errors arising 
from evaluation of £2^(2, A), Cz{z,X), Vg{z), and g {z) are all of size d^, and hence 
they too satisfy the required bounds. Similarly to ( |5.3| ), evaluation of and Dj\f 
yields errors of relative size S^, that is. 



(6.11a) 
(6.11b) 



computed Db 
computed 



Db- 



Gb, 



Gb = M^u, 
G^f = (5u/m, 



where Gb and G^ are diagonal matrices. 

We now obtain all the estimates in (|6.10 ) by a direct application of Corollary 4.4 , 
with the exception of the estimate for (As — As). Since the expressions for recovering 
As are identical in procedures condensed and augmented, we can apply expression 



(5.12) from Section 5.1 to deduce that the desired estimate holds for this component 
as well. □ 



The only difference between the error estimates of Theorem 5.1 for the condensed 



system and those obtained above for the augmented system is that the A\j\f com- 
ponents arc slightl y le ss accurate in the augme nted case. If we work through the 
analysis of Section 5^ with the estimate ( 6.10c ) replacing ( 5.13c ), we find that the 
main results are unaffected. Therefore, we conclude this section by stating without 
proof a result similar to Theorem 5.3. 
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-1.0 


-0.9 


-1.9 


-1.9 
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1 -2.7 


-1.5 


-1.3 


-1.2 


.9193 


(0.99, .19) 


5 -9.4 


-6.7 


-6.3 


-4.6 


1.0 


(1.04,. 23) 


6 -11.4 


-8.7 


-8.3 


-5.9 


1.0 


(1.04, .23) 


7 -13.4 


-10.7 


-10.3 


-3.8 


.9999 


(1.04,. 23) 


8 -15.4 


-12.7 


-12.3 


-1.2 


.9439 


(1.04,. 23) 


9 -17.1 


-13.9 


-13.4 


-0.6 


.9723 


(1.10, .20) 



Table 7.1 

Details of iteration sequence for PDIP applied to (2.8), with steps computed by solving the 
augmented system. 



Theorem 6.3. Suppose that all the assumptions of Theorem |5.q hold, except 
that the step (Az, AA, As) is calculated by using the procedure augmented with a 
factoriza tion/ solution algorithm that satisfies Condition 6A_- Then the conclusions of 
Theorem p.q hold. 

7. Numerical Illustration. We illustrate the results of Sections || and ^ using 
the two- variable example (2.8). Consider a simple algorithm that takes steps satisfying 
(3.8) with t set rather arbitrarily to t = /J,^e. (The search directions thus used are 
like those generated in the later stages of a practical primal-dual algorithm such as 
Mehrotra's algorithm ||l^.) We start this algorithm from the point 

zo = (1/30,1/9)^ Ao = (1,1/5)^, So = (1/10, 1/2). 



(It is easy to check that the conditions (3.11) are satisfied at this point for a modest 
value of C .) At each step we calculated dmax, defined in Section 5.3, and took an 
actual step of .99Q;max. 

We programmed the method in Matlab, using double-precision arithmetic. In our 
first experiment, we solved the formulation (4.1) of the linear equations using Matlab's 
standard Gaussian elimination solver for general systems of linear equations, which 
was analyzed in Section 6.4. From Theorem |6.2[ , the estimates (6.1C) apply to this 
case. 

Results are tabulated in Table 7.1. Note first the size of the component ||y-^AA/3||, 
which grows as /i decreases below u^/^, in accordance with ( ^.lOb ). (We cannot 
tabulate the difference ||y-^(AAB — AAgH) because of course we do not know the true 



step (Az, AA, As), but since the true step has size 0(/i) (Corollary 4.3), the error is 
dominated by the term y-^AAg in any case.) As predicted by (5.17), the maximum 
step dmax becomes significantly smaller than 1 as /i is decreased below u^/^. As 



indicated by ( 5.18D , however, good progress still can be made along this direction 
(in the sense of reducing /i and the norms of the residuals rf and r^) almost until ^ 
reaches the level of u. In fact, between iterations 5 and 8 we see the reduction factor 
of 100 that we would expect by moving a distance of .99 along a direction that is close 
to the pure Newton direction. The component with the large error — T^-^AAg — does 
not interfere significantly with rapid convergence, but only causes the A iterates to 
move tangentially to S\. This effect may be noted in the final iterate where the value 
of A changes significantly. In some cases, however, when the current A is near the 
edge of the set this error may result in a severe curtailment of the step length. 
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it6r 


log fl 


iog \\n.z\\ 


log 11 U ^ABll 


lOg\\V 1\Ab\\ 




A 





-1.0 


-0.9 


-1.9 


-1.9 


.9227 


(\ 00 20) 


1 


-2.7 


-1.5 


-1.3 


-1.2 


.9193 


(0.99, .19) 


5 


-9.4 


-6.7 


-6.3 


-4.6 


1.0 


(1.04,. 23) 


6 


-11.4 


-8.7 


-8.3 


-5.7 


1.0 


(1.04, .23) 


7 


-13.4 


-10.7 


-10.3 


-8.3 


1.0 


(1.04,. 23) 


8 


-15.4 


-12.7 


-12.4 


-10.3 


1.0 


(1.04,. 23) 


9 


-17.4 


-14.7 


-13.3 


-12.3 


1.0 


(1.04,. 23) 



lABLE I . Z 

Details of iteration sequence for PDIP applied to {2.S), with steps computed by solving the 
condensed system. 



Next, we performed the same experiment using the condensed formulation ( 3.10| ) 
of the linear system, as described in Section |^. Results are shown in Table 7.2. The 
main difference with Table 7.1 is that there is no increase in the value WV'^AXbW 



as /i approaches unit roundoff; this component appears to decrease at the same rate 
as the other step components. This observation can be explained by our analysis of 
the case in which the cancellation error term /g incurred in the evaluation of gsiz) 
satisfies (5.14). We calculated the product V'^{gB{z) + fs) (the product of V with 
our computed version of gi3{z)) and found it to be exactly zero on iterations 7, 8, and 
9. Therefore, using Taylor's theorem, ( 2.13| ), and Theorem we have 



V^fB = -V^gsiz) = -V^VgB{z*){z - z*) + 0(11^ - z*f) = O(fi^). 



Hence, ( 5.15| ) together with Corollary |4.3| shows that AXb — 0{fj.), which is con- 
sistent with the results in Table 7.2. Note too that because of the higher accuracy in 



the y^AAe component, the maximum step length stays very close to 1 during the 
last few iterations. By comparing Tables 7.1 and 7.2, however, we can verify that 



the convergence of /i to zero, and of the iterates to the solution set, is not materially 
affected by the presence or absence of the large error in AXb- 



To show that the lack of cancellation effects in Table 7.2 cannot be assumed in 
general, we modified problem (pM slightly, changing the second constraint to 



(7.1) 



(zi - V5y + z^ — — < 0. 



3\/5 



The primal and dual solutions remain unchanged, and we ran the condensed-equations 
version of the algorithm from the same starting point as above. Results are shown in 
Tabic 7.3. We observed that gsiz) did not escape cancellation errors in this instance 
and, as in Table 0, we observe significant errors in AXb that do not materially 
affect the convergence of the algorithm to the solution set. 

8. Summary and Conclusions. In this paper, we have analyzed the finite- 
precision implementation of a primal-dual interior point method whose convergence 
rate is theoretically superlinear. We have made the standard assumptions that appear 
in most analyses of local convergence of nonlinear programming algorithms and path- 
following algorithms, with one significant exception: The assumption of linearly inde- 
pendent active constraint gradients is replaced by the weaker Mangasarian-Fromovitz 
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1 T est* \ /~\ fV 1 1 

ITjcI lug fl 


log \ \l\Z\\ 


iog||ty L^AisW 


lr.(T WV'^ A\^\\ 
log II V AlAgl 




A 


-1.0 


-0.9 


-2.1 


-2.3 


.9161 


(\ 00 20) 


1 -2.7 


-1.5 


-1.3 


-1.4 


.8872 


(0.99, .20) 


5 -7.6 


-5.7 


-5.7 


-4.2 


.9999 


(.93, .29) 


6 -9.5 


-7.7 


-7.7 


-6.3 


1.0 


(.93, .29) 


7 -11.5 


-9.7 


-9.7 


-4.3 


.9999 


(.93, .29) 


8 -13.5 


-11.7 


-11.5 


-2.6 


.9960 


(.93, .29) 


9 -15.3 


-13.5 


-11.7 


-0.6 


.7386 


(.93, .29) 



Table 7.3 

Details of iteration sequence for PDIP applied to (2.i), (7.1), with steps computed from the 
condensed system. 



constraint qualification, which is equivalent to boundedness of the set of optimal 
Lagrange multipliers. Because of this assumption, it is possible that all reasonable 
formulations of the step equations — the linear system that needs to be solved to ob- 
tain the search direction — are ill conditioned, so it is not obvious that the numerical 
errors that occur when this system is solved in finite precision do not eventually ren- 
der the computed search direction useless. We show that although the error in the 
computed step may indeed become large as ^ decreases, most of the error is restricted 
to a subspace that does not matter, namely, the null space of the matrix V(7b(z*) of 
first derivatives of the active constraints. Although this error causes the computed 
iterates to "slip" in a tangential direction to the optimal Lagrange multiplier set, it 
does not interfere with rapid convergence of the iterates to the primal-dual solution 
set. 



We found that the centrality conditions ( 3.11 ) , which are usually applied in path 



following methods, play ed a crucial role in the analysis, since they enabled us to 



establish the estimates ( 3.1(: ) in Lemma 3^ concerning the sizes of the basic and 



nonbasic components of s and A near the solution set. The analysis of Section 



culminating in Corollary i.4 , finds bounds on the errors induced in step components 
by certain structured perturbations of the step equations. We show in the same 
section that the exact step is 0{fj,), allowing the local convergence analysis of Ralph 
and Wright p2[ to be extended from convex programs to nonlinear programs. 

In Sections ^ and ^ we apply the general results of Section ^ to the two most 
obvious ways of formulating and solving the step equations; namely, as a "condensed" 
system involving just the primal variables z, or as an "augmented" system involv- 
ing both z and the Lagrange multipliers A. In each case, the errors introduced in 
finite-precision implementation have the structure of the perturbations analyzed in 



Section ^ so the error bounds obtained in Corollary 4.4 apply. In Section 5.3 (whose 



analysis also applies to the computed solutions analyzed in Section |^) , we show that 
the potentially large error component discussed above does not interfere appreciably 
with the near-linear decrease of the quantities ^J■, rf, and rg to zero along the computed 
steps, indicating that until /i becomes quite close to u, the convergence behavior pre- 
dicted by the analysis of the "exact" algorithm will be observed in the finite-precision 
implementation. We conclude in Section [t] with a numerical illustration of our ma- 
jor observations on a simple problem with two variables and two constraints, first 
introduced in Section 3. 
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Appendix A. 



Justificat ion o f the Estimates ( 5.17 ) and ( 5.18 ). 

To pr ove ( 5.17 ), we use analysis similar to that o f S. Wright From the 

definition (3_^) of /i, and the centrality condition ( 3.11c ), we have that 

XiSi — Q{^), for alH = 1, 2, . . . , m. 



Hence, from the third block row of (3. J) and the assumption (3.7) on the size of i, we 
have that 



(A.1) ^+^^_i_4 

Ai Si SiXi 



-l + 0(/i). 



for alH = 1, 2, . . . , TO. 



We have from Lem ma |3.2| and (4.36) that AA^/A^ — 0(/i) for all i ^ B. Hence, by 
using ( |3.16a ) from (3.2) together with (A.l), we obtain 



(A.2) 



As, = 



for aU i e B. 



For the computed step components Asg, we have by combining (5.13a) with (A.2) 
that 



(A.3) 



As, 



-S,; + (5u + 0(Ai^), 



for aU i e B. 



Therefore, if Si + aAs^ = for some i B and some a G [0, 1], we have by using 
(3.16a) again that 



(A.4) 



(1 - a)s, = (5u + 0(^2) 

(1 — a) = Su/^J' + 0{fi), for any i £ B. 



Meanwhile, for i e Af, we have from Lemma 3.2, ( 4.36 ), and ( 5.13a ) that 

(A. 5) Si + aAsi > 0, for all i e J\f and all a e [0, 1], 

so the components AsAf do not place a limit on the step length bound dmax. For the 
components AAjv", we have by using Lemma 3.2, (4.36), (5.13c), and (A.l) that 

AAi = -Aj + fiSu + 0{fi'^), for all i e Af. 

There fore, if A^ + aAA,; = for some i e A/" and some a G [0, 1], we have by arguing 
as in ( |A^ that 



(A.6) 



1 - a (5u + 0(/i). 



Finally, for i G we have from Lemma 3.2 that A.; — 0(1), while from (4.36), (5.13a), 
and ( 5.13b| ), we have that 

(A.7) AA, = 0{fi), AA, = 0{fi) + d^/fi, for all i G B. 

Therefore, we have for u that 

(A.8) A, + aAA, > 0, for aU i G 6 and aU a G [0, 1]. 
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By combining the observations (A^), ( |A.5| ), ( [A. 61) , and (A^), we conclude that there 
is a value dmax satisfying 



e [0, 1], 1 - a„ 



such that 



(A, s) + a(AA, As) > 0, for aU a e [0 1 Q^inax I , 



proving the claim (5.17). By making various simplifications to the analysis above, it 
is easy to show that (5.19) holds as well. 

We now prove the claims ( ^.18 ) concerning the changes in the feasibility and 
duality measures along the computed step. 

From (1.2), ( [3.11a ), and the first block row of (3.8), we have 



(A.9) 



rf{z + aAz, A + aAA) 
^ Cz{z + aKz, A + aAA) 

= Cz{z, A) + aCz,Xz, \)Kz + aVg(z)AA + 0{a^\\Kzf) 
= (1 - a)£^(z. A) + aC,z{z, A)(A^ - Az) + aVg^{z){KXB - AAb) 
+aVgj^{z){K\^ ~ AAaa) + 0(a2||A^|P). 



From ( 4.3(: ) and ( 5.13a ), we have Az = (5u + 0(m), so for /i ^ u and a G [0, 1], we 
have 



(A.IO) 



a'||Az|p = 0(/i2 



From the definition (|2.13| ) of the SVD of Vgg(z*), Theorem U, and ( ^.13a]) , we have 
that 

V56(z)(AAf3 - AAe) = V.gg(z*)(AAs - AAg) + 0{\\z - z*||||AAg - AAe||) 
= C/EC/^(AAg - AAg) + 0{fi)6u/fi 
(A.ll) =<5u. 



Note that the larger error ( 5.13b| ) in the component ^-'"(AAe — AAe), which is present 
when MFCQ is satisfied but not when LICQ is satisfied, does not enter into the 
estimate (A.ll). By substituting this estimate into (A.9) together with estimates for 



Az — Az and AA^v' — AAa/' from (5.13), we obtain that 

rf{z + aAz, A + aAA) = (1 - a)r/ + + 0{fi^), 



verifying our claim (5.18c). The potentially large error ( 5.13b| ) does not affect rapid 
decrease of the rf component along the computed searc h direc tion. 

For the second feasibilit y mea sure , we have from ( 3.11b| ), the second block row 
of d"" 



and the estimates ( 5.13a ) and ( A.1C| ) that 

rg{z + aAz, s + aAs) 
= g{z + aAz) + s + aAs 

= g{z) + aVg{zfKz + s + aAs + 0(a^|lAz|p) 

= (1 - a)(.g(z) + s) + aV5(z)^(Az - Az) + a(As - As) + 0(^2) 

^{l-a)rg + 5^ + 0{y?), 
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verifying ( ^.18d ). 

To examine the change in fj,, we look at the change in each pairwise product A^s^ 
i — 1,2, . . . ,m. We have 

(\, + aA\^){s, + aAs,) 
= XiSi + a{siAXi + XiAsi) + a^AsiAXi 
(A.12) = Xis^ + a{siAX^ + X^As,) + as,(AAi - AA,) + aXi{Asi - Asi) 
-l-a^AAjAsj. 



From the last block row in (3.8), the estimate t = 0{^^) (^77|), and the estimate ( 4.36| ) 
of the exact step, we have 



(A.13) 



AiS, + a(siAA, + A,As,) = (1 - a)XiSi + 0{i?). 



From ( 4.36 ) and ( 5.13 ), we have 

(A.14) AA.ATs, = + 0{^i)){0{^i) + 6^) = 6^ + Oi^j.^), 

since /i ^ u. For i £ B, we have from Lemma ( 5.13a ), and ( 5.13b| ) that 

(A.15) Sj(AAj - AAi) = 0(a*)'5u/m = (^u, for alH e 6. 

For i € Af, we have from Lemma |3.2| and ( 5.13c| ) that 

(A.16) s,(AA, - AA,) = M<5u, foraUieW. 



For the remaining term Xi{Asi — As^), we have from Lemma 3.2 and (5.13a) that 
(A. 17) Ai(Asi - Asi) = (5u, for alH = 1, 2, . . . , m. 



By substituting (|A.13|) - (|A.17|) into (|A.12[) , we obtain 



(A.18|A, +aAA,)(s, + aAs^) = (1 - a)A,s, + + 0(^^), aU i = 1,2, . 



Therefore, by summing over i and using (^, we obtain ( |5l8b| ) 
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